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ABSTRACT 


This research demonstrates the Galerkin FEM's ability to provide approximate sol¬ 
utions of second order, nonlinear, one dimensional, two point boundary value problems. 
The research concentrates on the development of linearization, iteration, and interpo¬ 
lation strategies for the solution of differential equations containing the nonlinear u 2 
term. Additionally, various numerical considerations are explored. Over 2000 cases 
were analyzed using various strategies and results detailing the efficacy of strategy com¬ 
binations are presented. A linearization strategy known as quasilinearization consist¬ 
ently yielded excellent approximate solutions of the nonlinear differential equations 
investigated. It converged in a minimum number of iterations and was capable of solv¬ 
ing equations which have large function order and activity over their specified domain. 
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I. INTRODUCTION 


Differential equations which describe natural phenomena are developed through the 
application of the conservation principles of mass, momentum, and/or energy to a con¬ 
tinuous media. After application of these principles to a small or differential element 
of the media, the resulting equation provides a relationship between one or more deriv¬ 
atives of an unknown function whose behavior is desired. The differential equation is 
linear when it contains only derivatives of the function with either constant or inde¬ 
pendent variable function coefficients. The solutions of linear differential equations are 
unique based on their boundary and/or initial conditions. Various analytical and nu¬ 
merical strategies have been developed to obtain these solutions. The finite element 
method (FEM), which is used in this research, is one of the more popular and accurate 
techniques that has been developed to solve differential equations. 

When the coefficients of the derivatives are functions of the desired unknown, (that 
is, the dependent variable), or when the dependent variable and/or its derivatives do not 
appear linearly in the differential equation, then the differential equation becomes non¬ 
linear and conventional analytical techniques usually do not work. Nonlinear operators 
occur in a number of differential equations that describe the behavior of natural phe¬ 
nomenon, e.g., the Navier-Stokes equations for fluid flow, a beam on an inelastic foun¬ 
dation, the Falkner Skan equation, etc. This research investigates various strategies for 
solving nonlinear, second order, one dimensional, two point boundary value problems 
using FEM analysis. The Galerkin method of weighted residuals (MWR) with discrete 
basis functions is the FEM technique used to compute approximate solutions of these 
equations. 

A. LINEAR DIFFERENTIAL EQUATIONS 

Under most circumstances, this particular FEM provides excellent approximations 
to linear differential equations of the form 

2u-f= 0, 0 <x<D (1.1) 

with appropriate essential and natural boundary conditions where 

j ^2 

• if is the sum of arbitrary linear operators such as — (), (), etc. 

• x is the one-dimensional independent variable 
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• u is the dependent variable of x 

• / is an internal excitation to the system 

• D is the domain of the differential equation 

If the dominant operator in equation (1.1) is of odd order (non self-adjoint form), then 
the Galerkin method doesn't work and the Petrov-Galerkin FEM must be utilized. The 
Galerkin FEM transforms the differential equation into a system of linear algebraic 
equations of the form 

Au — F = 0 (1.2 .a) 

where 

• A is an NxN coefficient matrix characterizing the operator(s) if, and N is the 
number of system degrees of freedom (DOF) in the approximation. 

• u is the Nxl FEM approximate solution vector 

• F is an Nxl FEM vector representation of the excitation function/ 

This system of equations is readily solved by matrix algebra for the response variable, 

u = A _1 F (1.2.6) 

where equation (1.2.b) implies the solution of equation (1.2.a) and does not mean that 
A inverse is actually obtained. 

Chapter II of the research demonstrates the Galerkin Method's ability to accurately 
model linear differential equations, regardless of the magnitude of the solution, i.e., 
u — x 2 , x \ etc. Additionally, various methods of modeling the excitation function, /, are 
examined. 

B. NONLINEAR DIFFERENTIAL EQUATIONS 

The remaining research examines various nonlinear, second order, differential 
equations of the form 

ifw + if(u)-/= 0 0 <x< D (1.3) 

with appropriate boundary conditions where all terms are as previously defined and 
if (u) is the sum of arbitrary nonlinear terms such as u ~ , sin(u), u 1 , etc. Though this 
research concentrates on the u 2 nonlinear term, the principles presented should allow for 
the adequate analysis of any nonlinear term that might be encountered. 




The first step in the solution of equation (1.3) is linearization of the nonlinear term. 
Once the equation is linearized with respect to the dependent variable, the Galerkin 
FEM is used in an iterative fashion to solve a linear system of algebraic equations in the 
form 

Au + L*u — F = 0 (1.4) 

where the coefficients of L* are functions of u evaluated at the previous iteration, de¬ 
noted u,.„ and where subscripts i and i-1 refer to the present and past iterations. To 
begin the iteration process, a value is initially assumed for U;., and the system of 
equations is solved for u,. The new values of u ( and the input values of u M at each node 
are compared and tested against a convergence criteria. If convergence is not obtained, 
the new value of is substituted back into the system of equations for u M and the iter¬ 
ation process continues until convergence is obtained. The final iteration yields the 
FEM approximate solution of the nonlinear differential equation. 

Chapter III provides the general principles and considerations which are involved 
in solving nonlinear differential equations. Chapter IV utilizes these principles in the 
solution of two different equations containing the nonlinear term u 1 . Final conclusions 
are made based on the solution results as to which problem solving strategy yields the 
most accurate approximation while using the least amount of computer time. 
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II. LINEAR APPLICATIONS OF THE GALERKIN FINITE ELEMENT 

METHOD (FEM) 

This section examines the Galerkin method of weighted residuals (MWR) using 
discrete linear shape functions, that is the Galerkin FEM, as applied to linear second 
order, one dimensional, differential equations, that is, two point boundary value prob¬ 
lems. In section A, a general FEM procedure for the differential equation 

u" -f(x) = 0 a < x < b (2.1) 

with two boundary conditions, one at each end point of the domain, is presented. In 
particular, several strategies for the FEM representation of the excitation function f(x) 
are developed. 

In section B, the various strategies for FEM representation of the excitation func¬ 
tion are implemented on equation (2.1) for three different excitation functions. The 
three f(x) functions were selected to provide the solutions of x 2 , x\ x* to equation (2.1). 

A. GENERAL FEM PROCEDURE 

Here, a general procedure for the FEM formulation of equation (2.1) is presented. 
The Galerkin FEM process transforms a differential equation into a system of linear 
algebraic whose solution is an approximation to the exact solution of the differential 
equation. The transformation requires the following three steps: 

• Step 1: Form an N degree of freedom approximation, say u = G r u , where: 

u is the approximate solution 

G r is the lxN transpose vector of linear shape functions 

u is the Nxl vector of the FEM solution, u, at each node 

• Step 2: Form the Residual R-u”- f(x), where: 

f{x) is the excitation function in the differential equation 

• Step 3 Form the Galerkin integral equations J b QRdx = 0 

The evaluation of these integral equations gives a system of N linear algebraic 
equations whose solution is the approximation, u. The details of the FEM formulation 
for equation (2.1) follows. Substitution of the residual R into the Galerkin integral 
equation and separating terms yields 


4 





J*G Z"dx - J*G f{x)dx = 0 


(2.2.0) 


From step 1, u" = (G T u)" = (G T )''u. Substitution of this into the first term of equation 
(2.2.a) and moving the second term to the right hand side yields 


J fc G(G T )''<fr u = J^G f{x)dx 


(2.2 .b) 


Note that the term on the left hand side is the FEM representation of ifu = u" and the 
term on the right hand side is the FEM representation of /( x). The result of an inte¬ 
gration by parts of the left hand side leaves 


G(G t )'u | * - J*G'(G T )'dr u = j b Gf(x)dx 


(2.2.c) 


Each term is now evaluated individually. 


Boundary Term . G(G T )'u|* 

Differentiating the equation in step 1 of the formulation process yields (G T )'u = u'. 
Substitution of u' into the above equation and evaluating it at the upper and lower limits 


G,{b)Z\b) 

C 2 (b)Z'{b) 


Gj(a)u'(a) 

G 2 (a)u'(a) 


G«' | b ~ Gw' | a — 


(2.3 .a) 


G^{b)u\b) 


G N (a)u'(a) 


In equation (2.3.a), G , denotes that it is the linear shape function associated with the /* h 
system degree of freedom (SDOF) at the /* h system nodal point. The ' a ' and ' b ' ar¬ 
guments of G and i/ are the endpoint values of the domain where these functions are 
evaluated. These endpoin .^ could be denoted by their system node identities, which are 
1 and N, at the left and right end points respectively. With this notation, equation 
(2.3.a) becomes: 
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Gm'L-Gm'1. 


G x {N)u'{N) 


G,(1)S'(1) 

G 2 (N)a'(N) 

— 

G 2 (1)S'(1) 

g n (AO«'(M 


G n (1)«'(1) 


(2.3 -b) 


Due to the Kroenecker delta propertyi of the linear shape functions, all terms are equal 
to zero except G s {N)u’(N) and G,(l)u'(l). This yields a single vector comprised of the 
natural boundary conditions at each end of the domain and is designated B where 


Gu' | b — Gu' | a = B 



0 

Z‘(b) 


(2.3.c) 


When natural boundary conditions, u\ are present, they are used for the corresponding 
u' in equation (2.3.c). 


Differential Operator , /*G'(G T )’</jr u 

The integral f‘G'( < G [ )'dx , associated with the u" operator, is reduced to the element 
coordinate level for evaluation and becomes the 2x2 element a matrix 




(2.4.fl) 


where £ is the local coordinate axis, l, is the length of the element, and g is the 2x1 
vector of linear shape functions, given by 


l GU) 


{o ! 


=; 

*j 
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i-4- 


(2 M>) 


as shown in Figure 1. Differentiation of the g functions gives the 2x1 g' vector, 


Substituting g' and (g T )' into the integral of equation (2.4.a) gives 


(2.4.c) 


+ _L _ J_ 

,2 ,2 


_ ' 4 r ± i 6 

' + JL L 4 <JH 0 -± +-L 

J ° e 4 2 


-it 


+i -r 
-i +i 


(2-4.rf) 



■o-l) 


g! ’ T 


Figure 1. One Dimensional Linear Shape Functions 


After construction in an element DO loop, the 2x2 a matrix for each element is then 
distributed into the NxN system A matrix in accordance with a correspondence table, 
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which relates local DOF to system DOF, and where N is the number of system degrees 
of freedom. Each element has two local degrees of freedom, LDOF 1 at the left end of 
an element, and LDOF 2 at the right end of the element. The correspondence between 
LDOF 1 and LDOF 2 of the /'*’ element and the /" system degree of freedom (SDOF J) 
is J = (/ — 1) + k where k is 1 or 2. Upon assembly of all the element matrices, we obtain 
the system Au term. 


Excitation , f b G/(*) dx 

This term determines the forcing function (or excitation) vector F , that is 


F = 



(2.5) 


and is obtained by assembling the 2x1 element excitation vectors f. The f vectors can 
be either modeled as a lumped approximation term in several different ways or integrated 
exactly to yield a consistent forcing function. In this study, two lumped approximations 
and the exact integration are developed and thereafter compared to determine which 
yields the more accurate solution. A third approximation method is also described. 
Although this third method is not used in the evaluation of the excitation function in this 
chapter on linear systems, it is used in the next chapter on the nonlinear systems portion 
of the research. It should be noted that as the number of elements approaches infinity 
and the element length approaches zero, each approximation technique described below 
yields the exact value of the excitation integral. 

• Midpoint Lumped Approximation of f(x) 

The midpoint lumped approximation method for evaluation of the excitation 
vector is the simplest and crudest approximation. This approximation involves 
evaluating the function f{x) at the midpoint of the element and multiplying this 
value by the element length. Half of this area (/(/,/2)/,) is then placed at the left 
element node and the other half at the right element node as illustrated in 
Figure 2 on page 9 for three different arbitrary f(x). For the monotonically in¬ 
creasing function in Figure 2.a, this method places too much area at the left local 
nodal point (LNP) and not enough at the right. Conversely, when f(x) is 
monotonically decreasing, too little area is placed at the left node while the right 
node gets too much. When f(x) is concave over the element, too little area is 
placed at each node (Figure 2.b), while for the convex case (Figure 2.c), each LNP 
receives too much area. 

• 1/4 - 3/4 Lumped Approximation of f{x) 

This technique, which is a refinement of the previous one, evaluates/( x) at the 
quarter point and three quarter point of an element. Each value is multiplied by 
half of the element length and the resulting areas are placed at the left and right 
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clement nodal points respectively as shown in Figure 3 on page 10 for the same 
three arbitrary functions. This method provides a better approximation than the 
lumped midpoint technique, especially for the monotonically increasing function in 
Figure 3.a, as it uses two points to capture the behavior of the curve instead of one. 



Figure 2. Illustration of Midpoiut Lumped Approximation 


• Linear Approximation 

1 his method approximates f(x) over the element in terms of the linear shape 
functions where 

/to * /to -/to)^i - +/to)^ j~ ^ (2-6) 

as shown in Figure 4 on page 10. It overestimates the area for concave curves 
(Figure 4.a) and underestimates for convex curves as shown in Figure 4.b. Note 
that/(jr,) and /(or,) can be generalized to coefficients /, and^ to give a better linear 
lit of the curve. Also, the linear shape functions can be replaced by higher order 
shape functions which provides an order approximation of/( x) as opposed to 
a linear one. This approximation method is not utilized for the evaluation of 
forcing functions, but is used later in the nonlinear portion of the research to ap¬ 
proximate other types of functions. 

• Consistent 

The consistent solution requires transforming/(*) into the element coordinate 
system, /(£), by substituting or = <*,+ <*; into /( x) and performing the integration 
over the lengths of the elements. The coordinate transformation is illustrated in 
Figure 5 on page 11. 
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Figure 5. Global to Local Coordinate Transformation 

The integration for each of the above methods is performed on the clement level, 
producing the 2x1 f vector. The f vector for each element is then distributed into the 
Nxl system force vector F in accordance with the local to global nodal point coirc- 
spondcncc. 

Substitution of A, F, and B into equation (2.2.c) yields a linear system of algebraic 
equations in the form 

— Au = F — B (2.7.o) 

The (F — B) term can be combined into a single vector designated F m leaving 

- Au = F m (2.7. b) 

which can be solved for u, the FEM approximate solution at each nodal point. 

B. RESULTS 

In order to obtain specific results, the following equations are analyzed over the 
domain 0 < jc < 2: 


u” = 2; 

u(0) = 0, u'(2) = 4 

u exact x 

(2.8) 

t" = 6*; 

u(0) = 0, u'(2) = 12 

u =x 3 

u exact ■* 

(2.9) 

' = 12jt 2 ; 

t/(0) = 0, u'(2) = 32 

Uexact ~ x 

(2.10) 
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The excitation function in each of the above equations is evaluated using the three 
method* previously described in 11.A. The detailed formulation of these different f vec¬ 
tors is shown in Appendix A for the forcing function of equation (2.9) where f{x) = 6x. 
The other two f(x) are evaluated in a similar manner and the results for all three f(x) 
functions are shown in Table 1. 


Table 1. FORCING FUNCTIONS FOR VARIOUS F(X) 


/(*) 

Midpoint Approx¬ 
imation 

1/4 - 3/4 Approxi¬ 
mation 

Consistent 

2 


[l. 1 
A 



a 

A 



A 

A 


6x 

3/. 


I 

3/. 

BJ 

s 

■ 

r 3 ai.+in 

|_3a/, + 2/,j 

I2x 2 

61. 

I 

f * '9m A B 

1 

6/. 



i 

f 6a J /. + 4a/. 2 +/. 3 l 
[6 a 2 !. + 8a/, J + 3/. 5 J 


The FORTRAN programs and results for a ten element analysis of each equation 
with the various forcing functions are provided in Appendix B. The first problem con¬ 
sidered is that presented by equation (2.8). Due to the nature of the forcing function in 
this equation, i.e., a constant, all three formulation options provide the same result. The 
FEM approximation, shown in Figure 6, provided the exact solution at each nodal 
point for a 10 element analysis. 

The solutions of equations (2.9) and (2.10) were then considered. The different 
forcing function formulations in Table 1 were used in solving these differential 
equations. For clarity purposes, only a portion of the plots comparing each solution 
process to the exact solution in an area of rapidly changing gradient are shown in 
Figure 7 on page 14 (for u" = 6.x) and Figure 8 on page 15 (for u” — 12or : ). 

The midpoint approximation method for F provides a solution which is larger than 
the exact at each nodal point for both equations (2.9) and (2.10). The approximation 
is worse for equation (2.10) in Figure 8 as the midpoint method provides an overesti- 
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Figure 6. Comparison of Solutions for Equation 2.8 Using 10 Elements. 

mation of the actual value of the excitation integral at each node due to its inability to 
account for the quadratic nature of the forcing function. 

The quarter,'three-quarter point approximation ofF provides solutions for equations 
(2.9) and (2.10) that are quite close to the exact solution at each nodal point. This 
technique provides a much better approximation of the area under the forcing function 
curve because it discretizes the area into two independent sections, where as the mid¬ 
point technique did not. Thus, this technique is quite accurate in approximating 
excitation functions such as jc j which are monotonic and quadratic in nature. 

The consistent formulation method provides the exact solution for both equations 
(2.9) and (2.10) at each nodal point, even when the solution curve has a rapidly changing 
gradient such as u * x* as shown in Figure 8. It was expected that this technique would 
provide the most accurate solution for all cases as it yields the exact area given by the 
Galcrkin excitation integral. 
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Figure 7. Comparison of Solutions for Equation 2.9 Using 10 Elements. 
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C. CONCLUSIONS 

Linear shape functions provide an efficient and accurate interpolation for approxi¬ 
mating second order, linear, one dimensional, differential equations using the Galerkin 
FEM, regardless of the magnitude of the solution. Therefore, they are utilized 
throughout the remainder of the research when approximating linear and nonlinear op¬ 
erators of second order or less. An open question which remains is whether higher order 
elements will work for nonlinear two point boundary value problems when linear ele¬ 
ments do not. 

Additionally, two important observations regarding the use of a consistent forcing 
function analysis are made. 

• The use of this technique in evaluating the excitation integral provides for a very 
accurate solution process. Thus, this method is employed in the remainder of the 
research whenever possible. In those cases where the integration cannot be per¬ 
formed analytically, a Simpson's Rule approximation to the integral is used so that 
the resulting error is kept to a minimum. 

• For linear problems, this method provides very accurate approximations over large 
domains using a minimum number of elements. To illustrate this point, equation 
(2.10) was solved over the domain 0 < x < 10 using only two elements. The FEM 
approximation provided the "exact" solution at jr = 5 and jc = 10 as shown in 
Figure 9 on page 17. 
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Figure 9. Solution of Equation 2.10 Over a Large Domain Using Two Elements 
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III. ONE DIMENSIONAL SECOND ORDER NONLINEAR SYSTEMS 

This chapter outlines various solution strategies which are later analyzed on their 
ability to solve second order nonlinear differential equations of the form 

&u+ Jzf(u) — f(x) = 0 xeD (3.1.a) 

where S£ denotes linear operators and if() denotes nonlinear operators. Because 
equation (3.1.a) is nonlinear, a closed form analytical solution is, in general, not possible. 
Therefore numerical solutions are obtained by a variety of approximation techniques. 
This chapter sets c . a general procedure, consisting of three steps, for obtaining nu¬ 
merical solutions of equation (3.1.a). The three steps, when used with the Galerkin 
FEM, have the effect of transforming the original nonlinear differential equation into a 
system of linear algebraic equations. 

The first step in the procedure consists of a 'linearization' of the nonlinear : Sf(u) 
term(s). Thus, Jaf(u) is transformed to if'u where if' can be obtained by a number of 
different strategies, three of which are described in section 111. A. 

Once the nonlinear differential equation, equation (3.1.a), has been transformed to 
a linear differential equation of the form 

<£u + <£* u =/(*) (3.1.6) 

the second step is associated with how the if u term in equation (3.1.b) is evaluated in 
the FEM Galerkin integral equations. A number of interpolation procedures are de¬ 
scribed for this step in Section III.B. 

The third step in the solution procedure is defining the iterative process by which a 
solution of the linear algebraic equations developed by the Galerkin FEM is obtained. 
In particular, the efficacy of the iterative method involves two considerations. 

• the selection of an initial estimate to begin the iteration 

• the methodology for subsequent iterations 

Section III.C describes a number of iteration strategies. 

Thus sections III.A, III.B, and III.C cover the three basic steps in the solution 
procedure; that is, linearization, interpolation, and iteration. The selection of a partic¬ 
ular strategy within each of these steps defines a specific solution procedure which can 
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be utilized in approximating the solution of equation (3.1.a). In addition to these 
sections, section III.D discusses several numerical aspects which either affect, or are 
used in evaluating, the efficacy of a solution procedure. 

A. LINEARIZATION STRATEGIES 

The first step in analyzing nonlinear equations using the Galerkin FEM is 
linearization of the nonlinear term(s). Three different linearization strategies are inves¬ 
tigated in this research. There is no implication that these are the only strategies that 
exist. 

1. Constant Linearization 

The constant linearization method transforms the nonlinear term into a function 
of u' where 

l£{u) x <£*u = v(u ) (3.2 .a) 

As discussed in Chapter 1, the solution process for nonlinear equations using the 
Galerkin FEM is iterative in nature. On the first iteration, u' is set equal to an assumed 
value at each node and for subsequent iterations, each nodal value of u is based on the 
FEM approximation. Thus, v(u') is a known function evaluated at each node and is 
taken to the right hand side, leaving a linear equation of the form 

jS?u = f(x ) - v («*) (3.2 .b) 

As an example, consider the nonlinear term i/u. The constant linearization technique 
linearizes this term as v(u‘) = (u')V, where («')' is evaluated using finite difference 
techniques. 

2. Classical Linearization 

The classical linearization strategy linearizes the nonlinear term as a known 
function coefficient multiplying the dependent variable where 

i?(u) « SC'u = (m(u*)) u (3.3 .a) 

which in a sense maintains the functional nature of the dependent variable by allowing 
it to be kept on the left hand side of the differential equation in a linear fashion. As in 
the constant linearization strategy, m{u) is a known function coefficient where the values 
of u' are assumed for the first iteration and are based on the FEM approximation for 
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subsequent iterations. Substituting the results of equation (3.3.a) into equation (3.1.a), 
the linearized differential equation takes the form 

&U+ (m(u*)) u = f(x) (3.3 .b) 

Using the nonlinear term given as an example in the previous section, that is u'u, the 
classical linearization technique provides two alternatives. The first is (m(u’)) u — ( U 'Y U 
which keeps the full effect of the dependent variable, u . The second is (m(u')) u = u'u’ 
which maintains the full effect of the derivative term, 

3. Quasilinearization 

Quasilinearization is covered extensively by Bellman and Kalaba in [Ref. 1: p. 
36] where the nonlinear term is set equal to q[u) and linearized as 

&(u) * = q (u) + (u- u) q’(u) 

= (✓(«•)) «+(,(:•*)—V (u)) 

Comparing equation (3.4.a) with equations (3.2.a) and (3.3.a), it can be seen that 
quasilinearization is a combination of constant and classical linearization with the coef¬ 
ficient functions determined in a different manner. Substituting the results of equation 
(3.4.a) into equation (3.1.a), the linearized differential equation becomes 

Seu + (?'(u*)) u = f(x) - ( q (u ) - uq '(«*)) (3-4.^) 

As an example, consider the u 7 nonlinear term. The quasilinearization technique defines 
q(u) — u 3 , from which q(u) — (u) 1 and q’(u) = 2u’. Substitution of these functions into 
equation (3.4.a) yields a linearization of the form if’u = 2 u'u — (u') 2 . As u and u begin 
to approach the same value at each node during the iteration process, S£'u begins to 
approach the original nonlinear term, namely u\ 

B. INTERPOLATION STRATEGIES 

The second step of the solution process is the evaluation of the Galerkin integrals. 
Based on the type of linearization strategy utilized, one or both of the following integrals 
are obtained. 


Linearization Vector: 



(3.5.a) 
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Linearization Matrix: GG T h(u)dx u r,^ (3.5 .b) 

"£) 

The constant linearization strategy' results in a linearization vector where h(u) in 
equation (3.5.a) is replaced by v(u‘) from equation (3.2.a). The classical linearization 
technique yields a linearization matrix where h{u') in equation (3.5.b) is replaced by 
m(u‘) from equation (3.3.a). The quasilinearization method utilizes both integrals where 
h(u) becomes ( q(u *) — uq‘(u')) in equation (3.5.a) and q'{u) in equation (3.5.b). Since 
u ’ is derived from an FEM approximation utilizing linear shape functions, it varies line¬ 
arly between nodes over each element. Thus, when equations (3.5.a) and (3.5.b) are re¬ 
duced to the element level for evaluation, there are numerous interpolation strategies 
available to approximate h(u’). Here, a few strategies are discussed for each of the above 
integrals. 

1. Linearization Vector 

Reducing the integral in equation (3.5.a) to the element level yields 

[WVS (3-6) 

J o 

which is quite similar to the forcing function integral discussed in Chapter II. Three 
interpolation techniques similar to those used for the forcing function are examined in 
this research. The integral in equation (3.6) yields a 2x1 element vector denoted as f"; 
the superscript meaning that the vector changes with each iteration as the values of 
u are updated. The f' are then distributed into a system linearization vector denoted as 
F‘ in accordance with the local to global nodal point correspondence. 
a. Midpoint Lumped Approximation 

This method evaluates h(u ’) at the midpoint of the element and brings this 
value outside the integral as a constant. Since u’ varies linearly over the element, its 
value at the midpoint of the element is simply the average of the values of u at the two 
nodes of the element, (i£), and (r^Vi)< • Substituting this into equation (3.6) leaves 
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(“j)i + (“j+i)i \ 2 


where h () denotes that the function h is evaluated at the argument of ( ). 
b. 1/4 - 3/4 Lumped Approximation 

This method takes the h{u) term inside the shape function vector yielding 


r - 

Jo 


V"t 

h(u) ~ 


(3.8.c) 


In the first term, u is evaluated at IJA while in the second term it is evaluated at 3/./4. 
These values are again easily determined due to the linear variation of u over the ele¬ 
ment and are given by equations (3.8.b) and (3.8.c). 


“*( T ) “ ( u 7)> + \ ~ (“7)0 


4 ( u j )i + 4 ( U j+l)i 


(3.8 .b) 


“’() - (“j*)i + T (K+i)i ~ 


Substitution of these expressions for u into equation (3.8.a) gives 


(3.8.c) 
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where again, h () denotes that the function h is to be evaluated at ( ). 
c. Linear Approximation 

This method evaluates u in a linear manner over the entire element as a 
function of the element coordinate where 

«* = “ 4 * ) + ( u i+i)>( ) ( 3 -M 

Substitution of this expression into equation (3.6) yields 



(3.9.*) 


Since h{u') is no longer a constant but a function of <J, this integral must be evaluated 
for each specific /i(u'). Examples of each of these approximations is provided in Ap¬ 
pendix C. 

2. Linearization Matrix 

At the element level, the integral in equation (3.5.b) becomes 


fW h{u)di (3.10) 

J o 

which is similar to the Galerkin FEM differential operator integral discussed in Chapter 
II, with the gg T term producing a 2x2 element matrix. In order to preserve the 2x2 na¬ 
ture of this matrix, h(u) must remain as a single term multiplying each term in the ma¬ 
trix. Two techniques for evaluating this integral are examined in this research. The 
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resulting 2x2 element matrix is denoted as 1’ where the superscript again indicates that 
the matrix is changing with each iteration. These are then distributed into the system 
linearization matrix L’ in accordance with the local to global nodal point correspond¬ 
ence. 

a. Midpoint Lumped Approximation 

This method evaluates u at the midpoint of the element using the same 
process as in the linearization vector analysis and brings h(u’) outside the integral as a 
constant. Substitution of this expression into equation (3.10) gives 


(3.11) 


where h{) denotes that the function h is to be evaluated at ( ). 
b. Linear Approximation 

This method transforms u' into a linear function of the element coordinate 
as given by equation (3.9.a). This expression is substituted into equation (3.10) giving 

-—1 

L [ 1_ 7 7 ~7) + ( “ j+l)l (7))^ (3,12) 

<r 

m 

As for linear evaluation of the linearization vector, equation (3.12) must be evaluated for 
each h(u"). An example of each of these approximations is provided in Appendix D. 

C. ITERATION STRATEGIES 

Having evaluated the Galerkin integrals and developed a set of linear algebraic 
equations, the last step in the solution process is the determination of an iteration 
strategy. The main goal of an iteration strategy is to obtain a convergent approximation 
in a minimum number of steps while at the same time, keeping the computational effort 
of the iteration process to a minimum. Two of the factors which control the rate of 
convergence within a specific linearization strategy are 
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• how close the initial assumed values of u’ are to the actual solution 


• the nature in which the FEM solution , u, obtained at the end of each iteration, is 
utilized to obtain a value of u’ to begin each of the subsequent iterations 

1. Initial Iteration 

In order to begin any iteration process, an initial estimate, or guess, must be 
made as to the value of the variable which is to be determined. The relative accuracy 
of this guess with respect to the true solution of the differential equation greatly affects 
the ability of the solution process to converge as well as its rate of convergence. If the 
starting point for the iteration process is too far away from the actual solution, the 
likelihood for divergence or convergence to a nonsolution of the differential equation is 
very high. 

The first step in formulating an initial iteration strategy is developing an idea 
as to what the activity range of the solution of the differential equation might be, i.e., 
how much and how quickly is the dependent variable changing over the prescribed do¬ 
main. Since the solution is not known, this information must be obtained from the 
boundary conditions, the domain length, and insight as to the physics of the system. 

Two possible combinations of boundary conditions exist for two point bound¬ 
ary value problems. 

• Essential-Essential (E-E), where the magnitude of the dependent variable is speci¬ 
fied at each end of the domain 

• Essential-Natural (E-N), where the magnitude of the dependent variable is specified 
at one end of the domain, and the slope or rate of change of the dependent variable 
is specified at the other end. 

The essential-essential combination provides information as to the probable activity 
range of the solution over the system domain, which is referred to as function order in 
this research and is defined below. 

Function Order provides a relative magnitude of the range of values in the solution 
function as indicated by two essential boundary conditions. This relative 
magnitude is determined by writing the values of the essential boundary 
conditions in power ten exponent form and then taking the quotient of the 
maximum value over the minimum value. When the minimum valued 
boundary condition is 0.0, it should be written as 1 x 10° in order that the 
quotient does not become undefined. The magnitude of the exponent in this 
quotient defines the function order of the solution while the decimal value 
provides a ranking of how different function orders of the same magnitude 
compare. For example, take a differential equation which has boundary 
conditions of u(0) = 0 and u{ 2) = 25. These are written in pow’er ten expo¬ 
nent form as u(0) = 1.00 x 10° and u(2) = 0.25 x 10 2 . The quotient of these 
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two values is 0.25 jc 10 2 . Thus, the solution function is said to have a func¬ 
tion order of two. 

When an E-N boundary condition combination is specified, determination of 
the dependent variable activity range is not as straight forward as in the E-E situation. 
Instead, the activity range must be determined from a knowledge of the physics of the 
system as well as the actual values of the specified boundary conditions. The importance 
of properly estimating the activity range of the dependent variable cannot be overem¬ 
phasized as this estimation is utilized in determining an initial iteration strategy, which 
is the most critical step in the solution procedure. 

When the function order is zero or one, or the activity range is determined to 
be small based on the physics of the system, the magnitude of the dependent variable 
does not change appreciably over the prescribed domain. Thus, a reasonable estimate 
of the dependent variable for the first iteration would involve utilizing the essential 
boundary condition value(s). To examine the validity of this line of reasoning, the fol¬ 
lowing initial iteration strategies are utilized in the present research when the solution 
has a function order of one or less. 

• u* set equal to the value of the left essential boundary condition. 

• u set equal to the value of the right essential boundary condition 

• u set equal to the average value of the two essential boundary 7 conditions 

When the function order is greater than one, there is a chance that utilization 
of any of the above criteria for the initial iteration values could result in divergence or 
convergence to a nonsolution of the differential equation. Any number of guesses could 
be made for the initial iteration values, but the chances of a random guess providing a 
convergent solution of the differential equation are quite low. Instead, a systematic 
approach given by the following four steps is suggested. 

• Examine the physical system to which the differential equation applies and deter¬ 
mine which terms on the left hand side of the equation tend to dominate. 

• Neglect the nondominant term(s) and solve the equation for the dependent variable 
by any means possible, i.e., separation of variables, undetermined coefficients, 
successive integration, etc. 

• Determine the value of the dependent variable at each node and let these be the 
values used in the initial iteration. 

• If divergence results, reexamine the physical system and reevaluate the dominance 
of each of the terms. Then repeat steps two and three. 
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The process of determining initial iteration values that lead to a convergent solution can 
be very time consuming and frustrating. But, in most instances when the function order 
is greater than one, the selection of initial iteration values that are reasonably close' to 
the actual solution is crucial if a convergent solution is to be obtained. 

To amplify the above guidelines, consider a one dimensional constant cross- 
sectional area heat fin with nonlinear conduction given by 

k(r)r»-c 1 (7'-r oc ) = o (3.13) 

where c, is a constant based on the fin geometry and the convection coefficient, and T„ 
is the temperature of the convective fluid. When the fin is short or the value of k(T) for 
the given range of operating temperatures is much greater than c„ conduction tends to 
dominate and the convection term, c,(7 - T «,), is negligible. A value can then be as¬ 
sumed for k(7) based on the the physical system and boundary condition temperature(s). 
The equation thus becomes linear and can be integrated twice to yield T\x). T’ is then 
determined at each node which provides the values used for the initial iteration. 

When the constant c, is much larger than k(T), the convection term tends to 
dominate and conduction can be neglected. Thus, V is set equal to at each node for 
the initial iteration. In those cases where conduction dominates over one half of the 
domain and convection over the other, T' can be determined by using a combination of 
both these initial iteration strategies. The process of determining an initial iteration 
strategy is further developed in IV.A.5. 

2. Subsequent Iterations 

After the initial values of u are determined, the next step is to develop an iter¬ 
ation procedure for subsequent u values which will result in convergence with the least 
amount of computational effort. To aid in this development and subsequent discussions, 
the following notation is utilized. 

• («*), is the value of u used in the iteration process where j denotes the node number 
and i is the iteration number. 

• (u,), is the value of u returned by the FEM approximation where j and i are defined 
as above. 

Two different strategies for determining u are investigated in the present research. 
a. Previous Value Strategy 

This method takes the value of u from the previous iteration and sets it 
equal to u* for the next iteration process where 
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(3.14.0) 

This is the simplest iteration scheme, but does not take into account how u is changing 
during the iteration process. 

b. Average Value Strategy 

This method uses the average value of u from the last two iterations yielding 


, V + («,),—2 

(“i )i - 7 - 


(3.14.6) 


This method should enhance convergence when u is oscillating about the final conver¬ 
gent solution while at the same time not adversely affect those situations where u is 
converging monotonically. 

c. Additional Strategies 

Numerous other strategies can be utilized in an attempt to increase the rate 
of convergence. The following are not investigated in this research but are presented as 
topics requiring further research. 

• K-step Strategy - This method takes into account k previous iterations where 

« ( u y)i-l ( u j)i—2 "b ( u ])\-k x 

(Wj), -- 1 --- (3.14.C) 

• Weighted Average Strategy - A method which assigns a weighting factor to each 
iterative value of'u where 


. W,(«?j)i_, -H W 2 (S j ) i _ 2 

(“l 'i ~ w, + w. 


(3.14-0 


The objective is to find the optimum combination of weighting factors, w, and w 2 , 
that yields the minimum number of iterations. 

Weighted K-step Strategy - A combination of the previous two strategies where 

w l(“j).-l + w 2(“j),-2 + - + *]<(«]),-k 


(“j )i = ■ 


2 >. 


(3.14.e) 


Rate of Change Strategy - This method would require using some from of a Taylor 
series expansion to model how u changes from iteration to iteration. 


D. NUMERICAL CONSIDERATIONS 

The following numerical aspects are considered in evaluating the ability of each 
solution procedure to provide a convergent solution of the nonlinear differential 
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equations that are investigated. Again, a solution procedure is a specific combination 
of linearization, interpolation, and iteration strategies. 

1. Convergence/Divergence 

For a convergent solution, the absolute value of the difference between 
(Uj)t and (Uj),., at each node decreases as the number of iterations, i, increases. In order 
to achieve the best approximation for a given number of elements, it is desirable to let 
the maximum value of | (u^ — 1 at all nodes reach some minimum value before ex¬ 

iting the iteration process. Thus, the following percent difference convergence criterion 
is utilized: convergence is reached when the absolute value of the maximum difference 
between (U;)i and (w,)^ at any node divided by (uj, at the same node it less than .0001 
(.01%). This is shown mathematically in equation (3.15) and is but one of many con¬ 
vergence criteria that could be utilized. 


Convergence Criterion: 


MAX {(fyj-%.,1 

(«j)i 


< .0001 


(3.15) 


The solution procedure is considered a failure if any of the following situations occur. 

• convergence does not occur within 200 iterations 

• divergence occurs, i.e., it increases without bound as the number of iterations in 
creases 

• convergence to a nonsolution of the differential equation 


2. Critical Number of Degrees of Freedom (DOF) 

In some instances, especially when the function order is greater than two or the 
activity range is large, there may be a specific or critical number of DOF below which 
the solution procedure will not converge. T o examine this phenomenon, each solution 
strategy is initially evaluated using three DOF i.e., two elements. The number of DOF 
is then increased until either a grid independent solution is obtained or the number of 
DOF exceeds 100; and a critical number of DOF, if it exists, is determined. 

3. Stability 

Stability in numerical analysis applications is based on the approximation 
method's ability to converge as it relates to the time and displacement discretization that 
are utilized. For example, the finite difference explicit method is stable for r < where 
r = , while the Crank-Nicholsen implicit method is stable for any value of r, i.e., 

unconditionally stable. For this research, three types of stability are defined. 
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• Unconditional Stability - There is no critical number of DOF required to guarantee 
a convergent solution. 

• Conditional Stability - There is a critical number of DOF below which a convergent 
solution cannot be obtained. 

• Unstable - The solution process diverges regardless of the number of DOF utilized. 

4. Multiplicity of Solutions 

Nonlinear differential equations do not necessarily have a unique solution. The 
equations utilized in this research w r ere developed based on known solutions. When a 
solution different from them is obtained, that solution is checked at two different points 
in the domain by passing a parabola through the point in question and tw T o adjacent 
points. If the equations developed from both of these parabolas satisfy the differential 
equation, the solution procedure is considered to have provided a valid approximation. 

5. Boundary Condition Effects 

Two point boundary value problems require a boundary condition at each end. 
Two combinations are valid; either an essential (Dirichlet) boundary condition at each 
end or an essential at one end and a natural (Cauchy) at the other. The effect that each 
of these combinations has on the performance of each solution procedure is investigated. 

6. Computational Efficiency 

This research defines the most computationally efficient solution procedure as 
the one which provides the most accurate results for the least amount of computer time. 
CPU time by itself though is not an effective measure of efficiency. For example, a 
solution process that uses two elements will take much less time to run than one with 
40, but the 40 element solution is likely to be much more accurate. There should be 
some critical number of elements for a specific solution procedure beyond which the in¬ 
crease in CPU time for the additional calculations is not matched by a proportional in¬ 
crease in solution accuracy. 

The different solution strategies in this research are compared using a factor 
defined as 

CPU* = CPU x average % error (3.16.a) 

CPU is the amount of computer time (in seconds) required to complete the iteration 
process. A timing subroutine installed in the FORTRAN solution program starts when 
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the iteration process is entered and stops when it is exited. Average % error is the av¬ 
erage percent difference between the FEM approximation and the exact solution given 
by 


avg % 


error = 


N 

?, 1 % error 1 1 

i=i 


N 


^- _ number of system 
degrees of freedom 


(3.16.6) 


As the CPU time increases due to an increase in the number of elements, the 
percent error should decrease and approach zero. Therefore, the solution strategy 
yielding the minimum CPU' is defined as the most computationally efficient for a given 
equation and nonlinear operator. 
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IV. APPLICATIONS 


A. PRELIMINARIES 

1. Equations, Domains, and Boundary Conditions 

This chapter evaluates the performance of various combinations of the 
linearization, interpolation, and iteration strategies previously developed in solving sec¬ 
ond order, nonlinear, one dimensional, differential equations containing the u 1 nonlinear 
term. In order to investigate the effect of the many numerical considerations described 
in the preceding chapter, two nonlinear differential equations over three domains with 
appropriate essential and natural boundary conditions are solved. The differential 
equations considered were 

u" - u = 6 - 9x 4 u exact = 3* 2 (4.1) 

u" + U= 60* + 100* 6 Uexact - 10x 3 (4.2) 

Domains and Boundary Conditions 

• Domain 1: 0< * < 1 

Eqn. (4.1) u(0) *= 0, u(l) = 3 or u'(l) = 6 
Eqn. (4.2) u(0)*=0, u(l) - 10 or u'(l) = 30 

• Domain 2: 0 < * < 2 

Eqn. (4.1) u(0) = 0, u(2) = 12 or «'(2) = 12 
Eqn. (4.2) u(0) = 0, i/(2) = 80 or«'(2)= 120 

• Domain 3: 0 < * < 5 

Eqn. (4.2) u(0) = 0, u(5) = 75 or u'(5) - 30 
Eqn. (4.2) u(0) = 0, u(5) = 1250 or u'(5) = 750 

In all cases the left end point of the domain has an essential boundary condition. Either 
an essential or natural boundary condition is provided at the right end of the domain in 
order to investigate the effect of the two different boundary condition combinations, 
namely essential-essential (E-E) and essential-natural (E-N) as discussed in III.D.5. 

Equations (4.1) and (4.2) were developed by starting with a known solution, 
u, xfet , and working backwards to form a second order nonlinear differential equation. 
The equations were kept simple, i.e., only one linear and one nonlinear term, due to the 
number of solution procedures that required evaluation, as well as the many numerical 
considerations involved. Still, equations (4.1) and (4.2) are viable representations of two 
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engineering phenomena that are described by second order differential equations, such 
as 

• An axially loaded bar embedded in a nonlinear elastic medium as shown in 
Figure 10 on page 34. 

• A constant cross-sectional area heat fin with internal heat generation and nonlinear 
convection as shown in Figure 11 on page 35. 

2. Related Engineering Phenomena 
a. Bar Problem 

Consider the bar problem of Figure 10 subjected to distributed force p{x), 
embedded in a nonlinearly elastic media which exerts an opposing distributed force 
proportional to the square of the displacement, u. A free body analysis of a differential 
element yields 

F(x) + dF(x) + p{x)dx — F(jr) — u 1 (x)dx = 0 (4.3 .a) 

Cancelling the F(x) terms, taking the applied excitation p{x) to the right hand side of the 
equation and dividing by dx gives 

“7“ ~ “V) = “ P(*) (4.3.*) 


From solid mechanics, the following relations are known, 

F = a A 


o — cE 


t = 


du 

dx 


(4.4.a) 

(4.4.*) 

(4.4.c) 


where F is the axial force, o is the axial stress, A is the cross-sectional area, e is the 
strain, and E is Young's Modulus. Substitution of equations (4.4.c) and (4.4.b) into 
equation (4.4.a) provides a relation for the force as 

F - EAu' (4.4.d) 

Differentiating equation (4.4.d) with respect to x yields 

-T- = EAu" (4.4.<?) 
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Substitution of equation (4.4.c) into equation (4.3.b) leaves 


EA u" - u 2 = - p{x) (4.5) 

which is similar in form to equations (4.1) and (4.2), if the stiffness, EA, is set equal to 
unity. 




Figure 10. Axially Loaded Bar Embedded in a Nonlinearly Elastic Material 


b. Heat Fin Problem 

Now consider the heat fin problem of Figure 11 on page 35 where q(.r) is 
the volumetric heat gencration/unit length and the heat transfer coellicient, h, is a linear 
function of the temperature, T , that is h = a T. An energy balance on the fin differential 
element yields 

Qcond ({(x)dx (q c 0n d 4" dQcond) ^Qconv ® (4.6.<t) 

Cancelling the q tW terms, dividing by dx, and rearranging terms yields 


dUcond 

dx 


dtjconv 

dx 


= 


(4.6.h) 


The following relations are known from heat transfer principles. 






(4.7.a) 
(4.7 .h) 

\\ here A c is the cross-sectional area of the fin and is constant. A, is the surface area of 
the fin and is written as P„r where P is the perimeter of the heat fin. is the ambient 
temperature of the convective media, and k is the thermal conductivity coefficient. 
Substitution of equations (4.7.a) and (4.7.b) into equation (4.6.b) and performing the 
appropriate differentiation yields 

kA f r»-hP(7-7' 00 ) = q(x) 


1, A « * 

*\cond 

fifOHV ~ ^\j(7 7^) 


l.ctting h = a 7', = 0 , and dividing through by kA c yields 


7 ” —r~* 7 2 = — q(.v) 


(4.S./>) 


kA c kA c 

which is similar to equations (4.1) and (4.2) if coefficient (aP/kA t ) is set equal to unity. 





Figure 11. Heat Fin With Nonlinear Convection 
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3. Function Order and Domain Nondimensionalization 

The exact solutions on which equations (4.1) and (4.2) are based, and the three 
domains of these equations previously described, were judiciously chosen to provide in¬ 
formation as to the effect of function order on the performance of the solution proce¬ 
dures investigated in this analysis. The function order of each equation was determined 
over each domain based on the specified E-E boundary conditions and is provided in 
Table 2. 

In order to determine whether an advantage is gained by the 
nondimensionalization of a domain, the following investigation was undertaken for both 
equations over domains two and three with equation (4.1) over domain three provided 
as an example. In this case, the dimensional variable jc was replaced by the nondimen- 
sional independent variable rj = xj5. Since dtj = dx/5 and x = 5r), equation (4.1) was 
transformed to 

u" — 25u 2 = 150 — 140625 m 4 0<»?<1 

where differentiation now is with respect to tj. Analysis of each nondimensionalized 
equation established that the transfer of 'domain activity' to 'differential equation ac¬ 
tivity' did not result in any computational gain. 


Table 2. FUNCTION ORDER OF EQUATIONS (4.1) AND (4.2) 


Equation 

Domain 

Function Order 

(4.1) 

One 

1 

Two 

2 

Three 

2 

(4.2) 

One 

2 

Two 

2 

Three 

4 


4. General Solution Procedure 

Subsequent analyses using the solution procedures discussed next, show that the 
efficacy of any particular solution procedure depends strongly on the function order of 
the problem being solved. The general solution procedure consists of the three steps; 
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linearization, interpolation, and iteration; shown in Figure 12 on page 38. First, 
equations (4.1) and (4.2) are linearized using the linearization strategies developed in 
III.A. Secondly, an interpolation strategy developed in III.B is used with the Galerkin 
FEM and the linearized differential equation is transformed into a set of linear algebraic 
equations. Finally, these algebraic equations are solved iteratively using various com¬ 
binations of iteration strategies developed in III.C. The results of each solution proce¬ 
dure, where a solution procedure is a particular set of selected strategies, are then 
compared and evaluated based on their performance in approximating the solutions of 
equations (4.1) and (4.2). 

5. Initial Iteration Strategy 

The initial iteration strategy is the weak link in the overall solution procedure. 
For function order one problems, the initial iteration strategy is simply to utilize the es¬ 
sential boundary conditions as described in III.C.1. This is the strategy adopted for 
equation (4.1) over domain one, as well as equation (4.2) over the same domain because 
the latter problem becomes function order two only at the right boundary point of the 
domain. 

Another initial iteration strategy is utilized for equations (4.1) and (4.2) over 
domains two and three because their function order is greater than one. The physical 
systems to which equations (4.1) and (4.2) apply are generic in nature and do not provide 
the necessary information for determining which of the terms on the left hand side of the 
equation dominate the behavior of the system. However, since each of these equations 
is based on a known solution, the terms on the left side can be written as specific func¬ 
tions of jc. This information is utilized to devise an initial iteration strategy for each 
equation over domains two and three based on the indicated dominance of the linear and 
nonlinear terms. Normally, the solution of the nonlinear differential equation is un¬ 
known. In these cases, the physical system must be scrutinized and a determination 
made as to which processes and their corresponding operators dominate over each part 
of the domain. Having once assessed the dominance of each operator over various parts 
of the domain, the following strategies can be utilized to generate the initial iteration 
values. 

• Where the u " dominates, the initial iteration vector is obtained by solving 
u" */(jf). 

• Where the u 1 dominates, the initial iteration vector is obtained by solving u 2 =f(x ) 
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SOLUTION PROCEDURE 


NONLINEAR DIFFERENTIAL EQUATION 
& u + &(u) — f = 0 


£ 


LINEARIZATION 

STRATEGY 


LINEAR DIFFERENTIAL EQUATION 
<£*u — f — 0 


£ 


FEM WITH 
INTERPOLATION 
STRATEGY 


RECURSIVE LINEAR ALGEBRAIC EQUATIONS 
(A + A* t _i)u t - F* t _j - F = 0 


% 


INIT. ITER. STRAT. 
SUBS. ITER. STRAT. 
(i = ITER. CNTER) 

CONVERGED SOLN 
AFTER n ITER'S 


Figure 12. Solution Procedure for Nonlinear Differential Equations 
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In equation (4.1), u„. e , = 3x J which yields u" = 6 and u 2 = 9x\ These functions 
are plotted in Figure 13 on page 40 for the domain 0 < x < 2 and in Figure 14 on page 
41 for the domain 0<x<5. From Figure 13, the u” term is slightly more dominant 
than the u 2 term over the first half of the domain. Over the second half of the domain, 
the u 2 term is clearly dominant. Therefore, an appropriate initial iteration strategy is to 
utilize a procedure which neglects the u 2 term in determining the initial values of (u,') 0 
over 0 < x < 1 and neglects the u" term over 1 <, x < 2. The first part of the initial value 
procedure involves integrating equation (4.1) twice without the imposition of any 
boundary conditions. The second part requires taking the square root of the right hand 
side of equation (4.1). Thus, the initial iteration values used in solving equation (4.1) 
over domain two are determined using equation (4.9). 




0 < x < 1 
1 <x< 2 


(4.9) 


From Figure 14 it is apparent that the u 2 term dominates over a majority of the domain. 
Thus, the u” term is completely neglected and the initial iteration values are determined 
using equation (4.10). 

(u,WK~‘l (4.10) 


Nonlinear differential equation (4.2) was developed using u„«, - lOx 3 which 
yields u" = 60x and u 2 = lOOx*. These functions are plotted in Figure 15 on page 42 for 
a domain of 0 < x < 2 and in Figure 16 on page 43 over domain 0 < x < 5. Examination 
of Figure 15 shows the u" term slightly dominating the u 1 term over the first half of the 
domain and the u 2 term clearly dominating over the second half of the domain. Thus, 
the same procedure for determining initial iteration values as was utilized for equation 
(4.1) over domain two is employed and the results are given by equation (4.11). 

f 10* 3 + -^-* 8 0 <x < 1 

0 .-] , -^ ( 4 *(() 

^/eox+ioo* 6 l£x<2 
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Over the domain 0 < x < 5, Figure 16 clearly shows the domination of the u 2 term over 
a majority of the domain. Therefore the u" term is neglected and the initial iteration 
values are determined using equation (4.12). 

(«j ) 0 = -J6Qx + I00x 6 (4.12) 

The efficacy of the above initial iteration strategies was confirmed by the analyses which 
were undertaken. 

B. CONSTANT LINEARIZATION 
1. Problem Formulation 

The constant linearization strategy, described in III.A. 1, transforms the nonlin¬ 
ear term, u\ into a constant linear term as shown in equation (4.13) 

u 2 x&'u=(uf (4.13) 

where u' is determined as outlined in III.C. This process results in a linear differential 
equation of the form 

vT - ± (u*) 2 +/(*) x e D (4.14) 

where '+' is for equation (4.1),is for equation (4.2) and f{x) represents the excitation 
function in each equation. The Galerkin FEM formulation process outlined in Chapter 
II transforms equation (4.14) into 


G(G T )'u|*- f G'(G t )'^u = ± f G{ufdx + f Gf{x)dx (4.15) 

.I# j d J/j 

where a and b in the first term represent the value of x at the left and right boundary 
of the domain, respectively. The left hand side of equation (4.15) is similar to that of 
equation (2.2.c) and upon evaluation yields B - Au, where the vector B is only present 
when a natural boundary condition is specified. The integrals on the right hand side are 
now evaluated. 

Linearization Vector , jG(u'fdx 

This integral is evaluated as outlined in III.B.l where h(u') in equation (3.6) is 
set equal to (u 'f. The detailed formulation of the 2x1 f‘ element vectors for the three 
different interpolation strategies is given in Appendix C with the final results shown in 
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Table 3 on page 45. The 2x1 f’ vectors are then distributed into the system linearization 
vector, F’ , in accordance with the local to global nodal point correspondence. The F* 
vector continuously changes during the iteration process as the values of u’ are revised. 


Table 3. CONSTANT LINEARIZATION ELEMENT VECTORS 


Interpolation 

Strategy 


Midpoint 

Approximation 


1/4 - 3/4 
Approximation 


^ KUKm). J 


1 

2 

A 

2 


(- 3 - W .+1 «..><)’ 




K)? , 

(«;),(«;..), («•-,)? 


Linear 

Approximation 

/. 

4 

6 + 12 

Wm). . (tv,)? 




L 12 

6 + 4 J 



Excitation , jGf{x)dx 

Transforming x to the element coordinate system, £, the element integral for f 
becomes 


(4.16) 

where a, is the distance from the origin of x to the origin of i of the element for which 
f is being evaluated. Equation (4.16) is evaluated using the consistent technique de¬ 
scribed in 11.A, which is detailed in Appendix A, where/() is replaced by the appropri¬ 
ate excitation function in equations (4.1) and (4.2). The resulting 2x1 element excitation 
vectors are presented in Table 4 on page 46. The f vectors remain steady (constant) 
during the iteration process and are distributed into an Nxl system force vector, F, in 
accordance with the local to global DOF correspondence. 
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Table 4. ELEMENT FORCE VECTORS FOR EQUATIONS (4.1) AND (4.2) 


/to 

f 

6 — 9.x 4 

mm 


■ 

■ 

3 

B 

a: 

]-f“ f 

K 

]~ 30 '* [J] 

60jc + 100 .x 4 

30»«/.[] 
100 a 3 //| 

] + 
a 

+ 50a 2 // 

+ 50a 4 / 

[!]- 

,[j] + lOOa 5 //^] + 125a 4 /, 3 [ 3 ] + 

100 , 25 

T~ a H6j + l6'H7j 


Final FEM Equations 

Substitution of the matrix and vector equivalents for each integral into equation 
(4.15) yields a system of equations given by 

B —Au = F‘ + F (4.17.a) 

where the ± sign in equation (4.15) is incorporated into F\ The natural boundary con¬ 
dition vector B, when present for an E-N problem does not change during the iteration 
process and is taken to the right side and subtracted from F yielding F m , where the m 
subscript indicates that the excitation vector is modified for the given natural boundary 
condition. The system of equations then takes the final form of 

—Au,= F*_j + F ra (4.17.*) 

where subscripts i and i — 1 refer to the iteration counter. Equation 4.17.b is solved it¬ 
eratively for u, with F* changing after each iteration, until convergence is obtained. 

2. Results 

a. General 

Equations (4.1) and (4.2) were each solved over domain one using 24 dif¬ 
ferent solution procedures while 12 different procedures were utilized for both domains 
two and three. The FORTRAN programs utilized for the constant linearization strategy 
are contained in Appendix E. A summary of the strategies utilized in each solution 
procedure is shown in Table 6 on page 49 and Table 7 on page 50 for equation (4.1), 
and Table 8 on page 51 and Table 9 on page 52 for equation (4.2), with the following 
performance information provided in the results portion of each table. 
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• Convergence (III.D.l) 

• Stability (III.D.3) 

• Number of iterations required to obtain convergence 

• Average percent error (III.D.6) 

• CPU' (III.D.6) 

The results for each solution procedure over each of the domains were obtained using 
the number of elements shown in Table 5. Though these number of elements are not 
necessarily the number required to obtain a grid independent solution for each solution 
procedure, they do provide a common baseline upon which the performance of individ¬ 
ual solution procedures for a particular problem may be compared. This is also the 
number of elements utilized in the table of results for each of the solution procedures 
developed by the other two linearization strategies. 


Table 5. NUMBER OF ELEMENTS UPON WHICH SOLUTION PROCEDURE 
RESULTS ARE BASED 


Equation 

Domain 

Number of Elements 

4.1 

1 

10 

2 

20 

3 

25 

4.2 

1 

20 

2 

40 

3 

50 


Two general observations can be made based on the results of Table 6 
through Table 9. One is that the constant linearization technique begins to fail as the 
function order over the given domain begins to approach two. This is clearly shown in 
Table 8 as only half of the solution procedures provided convergent solutions of 
equation (4.2), whose function order is just barely two over the domain 0 < x < 1. As 
soon as the domain length was increased to 1.1, these remaining methods diverged as 
shown in Table 9. This is most likely due to equation (4.2) becoming fully order two 
as the domain is increased past x = 1.0. In order to determine where the solutions begin 
to breakdown for equation (4.1), the domain was extended in 0.1 increments until all 
solution procedures diverged. The results of this investigation showed that some sol¬ 
ution procedures for equation (4.1) were able to provide convergent solutions over a 
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domain of 0 < jc < 1.7 , which approaches the region where equation (4.1) shifts in 
function order from one to two. This failure to provide convergent solutions of 
equations (4.1) and (4.2) is not due to a lack of elements used in the approximation, but 
is a characteristic of this particular linearization strategy. 

The second observation is that the number of iterations required to obtain 
convergence is strictly a function of the iteration strategies utilized, both initial and 
subsequent, and is independent of the interpolation strategy. Also, it was found that the 
number of iterations is independent of the number of degrees of freedom utilized. That 
is, the same number of iterations were required to reach convergence whether two or 20 
elements were utilized. The accuracy of the approximation, on the other hand, seems 
to be independent of the iteration strategy, and only a function of the interpolation 
strategy and the number of degrees of freedom. The effect of each problem parameter 
in the different solution procedures for both equations (4.1) and (4.2) over domain one 
is now examined, as no convergent solutions were obtained for domains two and three. 

b. Boundary Conditions 

The use of an essential boundary condition at both ends of the domain in¬ 
creased the rate of convergence by a factor of two to five over those strategies which 
utilized an essential and natural boundary condition combination. In fact, only three 
of the twelve strategies which used an E-N boundary condition combination provided 
convergent approximations of the differential equation. Suffice it to say, the essential- 
natural boundary condition combination does not produce very efficient results when 
utilized with the constant linearization strategy and therefore this strategy should not 
be used for E-N problems. That being the case, the remaining comments are directed 
at those strategies which utilized an essential-essential boundary condition combination. 

c. Initial Iteration Strategy 

Each of the initial iteration strategies, described in III.C.l, provided nearly 
he same results within a specific subsequent iteration and interpolation strategy com¬ 
bination when the function order was small, as in equation (4.1). But, when the function 
order begins to approach two, as in equation (4.2), the use of the right essential bound¬ 
ary condition for the initial iteration values led to divergence, while the other two strat¬ 
egies provided similar convergent results. 

d. Subsequent Iteration Strategy 

The use of the previous value strategy (IIl.C.2.a) generally led to conver¬ 
gence using less iterations than the average value strategy (lII.C.2.b). The main reason 
for this is most likely a result of the constant linearization strategy providing for 
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Table 6. SOLUTION PROCEDURES AND RESULTS USING CONSTANT 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAIN 
ONE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

K)o 


Stab. 

# of 

% 

CPU* 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 




E 

E 

H 

C 

S 

7 

0.48 

0.0112 

0 < x < 1 

Prev. 

Midpt. 

E 

E 

M 

C 

s 

7 

0.48 

0.0081 

Value 

E 

E 

R 

c 

s 

8 

0.48 

0.0129 




E 

X 

n 

c‘ 

.... 

200 - 

.... 





E 

E 

n 

c 

s 

7 

0.12 

0.0028 

0 < ;r < 1 

Prev. 

1 ,'4-3/4 

E 

E 

M 

c 

s 

7 

0.12 

0.0028 

Value 

E 

E 

R 

c 

s 

8 

0.12 

0.0032 




E 

X 

n 

c* 

.... 

200 - 

.... 





E 

E 

m 

c 

s 

7 

0.16 

0.0027 

0 < jc < 1 

Prev. 

Linear 

E 

E 

M 

c 

s 

7 

0.16 

0.0043 

Value 

E 

E 

R 

c 

s 

8 

0.16 

0.0043 




E 

X 

n 

c* 

—— 

200 - 

.... 





E 

E 

n 

c 

s 

9 

0.48 

0.0127 

0 <jc< 1 

Avg. 

Midpt. 

E 

E 

M 

c 

s 

KM 

0.50 

0.0079 

Value 

E 

E 

R 

c 

s 

9 

0.50 

0.0132 




E 

X 

m 

c 

s 

38 

0.27 

0.0270 




E 

E 

n 

c 

s 

9 

0.12 

0.0025 

0 < or < 1 

Avg. 

1/4-3/4 

E 

E 

M 

c 

s 

7 

0.12 

0.0022 

Value 

E 

E 

R 

c 

s 

9 

0.11 

0.0037 




E 

X 

a 

c 

s 

38 

0.31 

0.0286 




E 

E 

n 

c 

s 

9 

0.16 

0.0041 

0 < jc < 1 

Avg. 

Linear 

E 

E 

M 

c 

s 

7 

0.15 

0.0035 

Value 

E 

E 

R 

c 

s 

9 

0.17 

0.0046 




E 

X 

a 

c 

s 

38 

0.47 

0.0471 

LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R= right essential boundary condition; 

M * average of L and R; C = convergence; D = divergence S — unconditional stability; 
CS = conditional stability; U = unstable 

NOTES: , 

a - Oscillates between two nonsolutions of the differential equation 







































































































































































































Table 7. SOLUTION PROCEDURES AND RESULTS USING CONSTANT 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAINS 
TWO AND THREE 


Results 


Domain 

Solution Procedure 

Iter. 

Interp. 

B. < 

Z.'s 

(4)« 


Strat. 

Strat. 

Lt 

Rt 

0 <x< 2 

Prev. 

Midpt. 

E 

E 

(1) 

Value 

E 

N 

(1) 

0 < x<2 

Prev. 

1,'4-3/4 

E 

E 

(1) 

Value 

E 

N 

(1) 

0 <x<2 

Prev. 

Linear 

E 

E 

(1) 

Value 

E 

N 

(1) 

0 <x< 2 

Avg. 

Midpt. 

E 

E 

(1) 

Value 

E 

N 

(I) 

0 < x < 2 

Avg. 

1/4-3/4 

E 

E 

(1) 

Value 

E 

N 

(1) 

0 <x< 2 

Avg. 

Linear 

E 

E 

(1) 

Value 

E 

N 

(1) 

0 < x < 5 

Prev. 

Midpt. 

E 

E 

(2) 

Value 

E 

N 

(2) 

0 < x < 5 

Prev. 

1/4-3/4 

E 

E 

(2) 

Value 

E 

N 

(2) 

0 < x < 5 

Prev. 

Linear 

E 

E 

(2) 

Value 

E 

N 

(2) 

0<x<5 

Avg. 

Midpt. 

E 

E 

(2) 

Value 

E 

N 

(2) 

0 < X < 5 

Avg. 

1 ,'4-3/4 

E 

E 

(2) 

Value 

E 

N 

(2) 

0<x< 5 

Avg. 

Linear 

E 

E 

(2) 

Value 

E 

N 

(2) 



LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R * right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS = conditional stability; U = unstable 


NOTES: (1) - See equation (4.9) (2) - See equation (4.10) 
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Table 8. SOLUTION PROCEDURES AND RESULTS USING CONSTANT 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAIN 
ONE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

(«;)» 


Stab. 

# of 

% 

CPU’ 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 




E 

E 

n 

C 

S 

12 

12.9 

1.4595 

0 < x < 1 

Prev. 

Midpt. 

E 

E 

M 

C 

s 

13 

13.0 

2.0773 

Value 

E 

E 

R 

D 

mm 


.... 





E 

N 

a 

D 

u 


.... 





E 

E 

a 

C 

s 

12 

2.12 

0.2198 

0 < x < 1 

Prev. 

1/4-3/4 

E 

E 

M 

C 

s 

13 

1.99 

0.3380 

Value 

E 

E 

R 

D 

U 


.... 





E 

N 

a 

D 

u 


.... 





E 

E 

a 

C 

s 

12 

4.62 

0.4922 

0 < x < I 

Prev. 

Linear 

E 

E 

M 

c 

s 

13 

4.77 

0.6664 

Value 

E 

E 

R 

D 

mm 


.... 





E 

N 

n 

D 

u 


—— 





E 

E 

a 

C 

s 

17 

12.8 

1.7412 

0 <x< J 

Avg. 

Midpt. 

E 

E 

M 

C 

s 

18 

13.1 

2.8829 

Value 

E 

E 

R 

D 

u 


—— 





E 

N 

n 

msm 

s 

32 

.... 





E 

WM 

a 

c 

s 

17 

2.22 

0.3619 

0 <x< 1 

Avg. 

1/4-3/4 

E 

E 

M 

c 

s 

18 

1.97 

0.3799 

Value 

E 

E 

R 

D 

U 


.... 





E 

N 

n 

C* 

s 

32 

.... 





E 

E 

a 

C 

s 

17 

4.58 

0.7091 

0 < x < 1 

Avg. 

Linear 

E 

E 

M 

C 

s 

18 

4.89 

1.0410 

Value 

E 

E 

R 

D 

u 


---- 





E 

N 

a 

C* 

s 

32 

—— 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R= right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS = conditional stability; U = unstable 

NOTES: . 

a - Converges to a nonsolution of the differential equation 
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Table 9. SOLUTION PROCEDURES ND RESULTS USING CONSTANT 
LINEARIZATION TO SOLVL EQUATION (4.2) OVER DOMAINS 
TWO AND THREE 


Solution Procedure 


Results 


Domain 

Iter. 

Interp. 

B. ( 

Z.'s 

«)o 


Stab. 


Strat. 

Strat. 

Lt 

Rt 


Prev. 

Midpt. 

E 

E 

( 1 ) 

D 

mm 

0 < or < 2 

Value 

E 

N 

(I) 

D 

u 

0 <x <2 

Prev. 

1 ,'4-3/4 

E 

E 

( 1 ) 

D 

mm 

Value 

E 

N 

( 1 ) 

D 

u 

0 < x < 2 

Prev. 

Linear 

E 

E 

( 1 ) 

D 

um 

Value 

E 

\ 

( 1 ) 

D 

u 

0 <x< 2 

Avg. 

Midpt. 

E 

E 

( 1 ) 

D 

u 

Value 

E 

N 

( 1 ) 

D 

u 

0 < x < 2 

Avg. 

1/4-3/4 

E 

E 

( 1 ) 

D 

mm 

Value 

E 

N 

mi 

D 

mm 

0 < x < 2 

Avg. 


E 

E 

0 ) 

D 

mm 

Value 

Linear 

E 

N 

(i) 

D 

U 

0<x<5 

Prev. 

Midpt. 

E 

E 

( 2 ) 

D 

mm 

Value 

E 

N 

( 2 ) 

D 

mm 

0 < x < 5 

Prev. 

1/4-3/4 

E 

E 

( 2 ) 

D 

u 

Value 

E 

N 

( 2 ) 

D 

u 

0<x<5 

Prev. 

Linear 

E 

E 

( 2 ) 

D 

mm 

Value 

E 

N 

( 2 ) 

D 

mm 

0<x< 5 

Avg. 

Midpt. 

E 

E 

( 2 ) 

D 

U 

Value 

E 

N 

( 2 ) 

D 

U 

0 <x< 5 

Avg. 

1 IA 'll A 

E 

E 

( 2 ) 

D 

mm 

Value 


E 

N 

( 2 ) 

D 

u 

0<x<5 

Avg. 

Linear 

E 

E 

( 2 ) 

D 

mm 

Value 

E 

N 

( 2 ) 

D 

mm 


LEGEND: E * essential boundary condition; N = natural boundary condition; 

L * left essential boundary condition; R = right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS * conditional stability; U = unstable 


NOTES: (1) - See equation (4.11) (2) - See equation (4.12) 



























































































































































monotonic convergence when used with an E-E combination. The essential boundary 
conditions provide a range that the values of u should be between and thus provide 
strong guidance during the iteration process. The one exception to this involves those 
cases where a natural boundary condition is given at the right end of the domain. In 
these instances, only the average value iteration strategy resulted in convergent sol¬ 
utions. This E-N combination does not fix the value of the dependent variable at the 
right end, so the approximate solution most likely oscillates about the exact solution 
during the iteration process. The average value strategy tends to enhance the conver¬ 
gence of this type of approximation, which might explain why this strategy resulted in 
convergence of the E-N problem while the previous value strategy did not. 
e. Interpolation Strategy 

The 1/4-3/4 interpolation method consistently provided the most accurate 
approximations as indicated by the low average percent error values. This is because the 
solution function for both equations (4.1) and (4.2) is a polynomial function of x and the 
1/4-3/4 method was shown in Chapter II to provide better approximations of the inte¬ 
grals of these functions than the midpoint or linear techniques. 
f Overall Performance 

The overall performance, which factors both computational effort (i.e., CPU 
time) and solution accuracy, of a particular solution procedure is indicated by its re¬ 
spective value of CPU'; where the lower the value, the more efficient the solution pro¬ 
cedure. In Table 6, the CPU times upon which the CPU* values are based were all at 
or below the clock subroutine accuracy of ± .03 seconds. Thus, the CPU* values in this 
table should be used with caution. A comparison of the solution procedures using av¬ 
erage percent error values (% Dif in Table 6) indicates the previous value iteration 
strategy combined with the 1/4-3/4 interpolation strategy provides the most accurate 
solution of equation (4.1) over domain one. 

Due to the increase in the number of elements and iterations required for 
convergence of equation (4.2), more CPU time was required by the solution procedures 
in Table 8. Hence, the CPU* times in Table 8 are all based on CPU times much greater 
than the clock subroutine accuracy and therefore are all valid. Again, the combination 
of previous value iteration and 1/4-3/4 interpolation strategies provided for the most ef¬ 
ficient procedure. 

3. Conclusions 

The constant linearization strategy can generally provide approximations of 
nonlinear differential equations when the function activity is less than order two over the 
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given domain. The relative 'crudeness' of this linearization technique requires that as 
much information as possit concerning the actual value of the dependent variable over 
the given domain be known. Thus, the knowledge of an essential boundary condition 
at each end of the domain is almost a prerequisite to obtain a convergent solution of a 
nonlinear differential equation when using this linearization method. If these conditions 
are met, an iteration strategy utilizing the magnitude of the smallest valued boundary 
condition for the initial iteration and the previous value strategy for subsequent iter¬ 
ations should result in convergence with a minimum number of iterations and expend 
the least amount of CPU time for a given number of degrees of freedom. The 1/4-3'4 
interpolation strategy should yield the most accurate approximation as most solutions 
of engineering problems are monotonically increasing or decreasing functions, or at 
worst, convex or concave over the given domain. 

C. CLASSICAL LINEARIZATION 
1. Problem Formulation 

The classical linearization strategy transforms the u 2 nonlinear term into a linear 
term as described in III.A.2 and shown in equation (4.18) 

u 2 x&*u = u*u (4.18) 

w'here u* is determined as outlined in III.C. Substitution of equation (4.18) into 
equations (4.1) and (4.2) yields a linear differential equation of the form 

u"±u*u—f(x) xeD (4.19) 

where the '—' is for equation (4.1), the '+' is for equation (4.2) and f(x) is again the 
respective excitation function in each equation. The Galerkin FEM formulation process 
transforms equation (4.19) into 

G(G T )'u | * - f G'(G t )Vjc u ± f GG J udx u = f Gf(x)dx (4.20) 

a j d J[) Jd 

The first two terms on the left side of equation (4.20) again provide B — Au where the 
B vector is present only when a natural boundary condition is provided. The integral 
on the right side of equation (4.20), which was evaluated in IV.B.l, gives the system 
excitation vector, F. The only term remaining to be evaluated is the third term on the 
left side of equation (4.20), the linearization matrix integral. 
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Linearization Matrix . f GG'u'dx u 

This integral is evaluated as outlined in I1I.B.2 where h{u) in equation'*3.10) is 
replaced by u\ The detailed formulation of the element linearization matrices, 1, for the 
two different interpolation strategies is given in Appendix D, with the final results pro¬ 
vided in Table 10. The 1' are then distributed into the system linearization matrix, L", 
based on the local to global nodal point correspondence. The L’ matrix is updated 
during each iteration as the values of u’ are revised. 


Table 10. CLASSICAL LINEARIZATION ELEMENT MATRICES 


Interpolation 

Strategy 

r 

Midpoint 

Approximation 

H9 


■ 

Linear 

Approximation 

i, r 3(u*),+(u,\,)/ K), + K- 1 ), 1 

12 [ K), + Km), K), + 3Km), J 


Substitution of the matrix and vector equivalents for each integral into equation 
(4.20) yields a system of equations given by 

B-Au + L*u=F (4.21.a) 

where the ± sign in equation (4.20) is incorporated into L\ Again, F and B are com¬ 
bined to yield F m and the system of linear algebraic equations takes the final form 

(-A+On-F, (4.21.6) 

Equation (4.2l.b) is solved iteratively for u,, with L’ M being calculated using u’., values, 
until convergence is obtained. 

2. Results 

a. General 

Sixteen different solution procedures were utilized to solve equations (4.1) 
and (4.2) over domain one and eight procedures were utilized for each equation over 
both domains two and three. The FORTRAN programs for constant linearization are 
contained in Appendix F. The different strategies utilized in each procedure and the 
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corresponding results are provided in Table 11 on page 56, and Table 12 on page 57 for 
equation (4.1); and Table 13 on page 58, and Table 14 on page 59 for equation (4.2). 
The number of elements utilized in each solution procedure is shown in Table 5 on page 
47. 


Table 11. SOLUTION PROCEDURES AND RESULTS USING CLASSICAL 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAIN 
ONE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

K)o 


Stab. 

# of 

% 

CPU' 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 




E 

E 

n 

C 

S 

6 

0.24 

0.0048 

0 < x < 1 

Prev. 

Midpt. 

E 

E 

M 

C 

s 

4 

0.24 

0.0040 

Value 

E 

E 

R 

c 

s 

5 

0.24 

0.0032 




E 

N 

a 

c 

s 

12 

0.10 

0.0030 




E 

E 

n 

c 

s 

6 

0.16 

0.0032 

0 <x< 1 

Prev. 

Linear 

E 

E 

M 

c 

s 

4 

0.16 

0.0016 

Value 

E 

E 

R 

c 

s 

5 

0.16 

0.0038 




E 

N 

a 

c 

s 

12 

0.38 

0.0115 




E 

E 

n 

c 

s 

8 

0.24 

0.0048 

0 < x < 1 

Avg. 

Midpt. 

E 

E 

M 

c 

s 

5 

0.24 

0.0048 

Value 

E 

E 

R 

c 

s 

8 

0.24 

0.0064 




E 

N 

n 

c 

s 

12 

0.08 

0.0019 




E 

E 

n 

c 

s 

8 

0.16 

0.0027 

0 <x< 1 

Avg. 

Linear 

E 

E 

M 

c 

s 

5 

0.16 

0.0043 

Value 

E 

E 

R 

c 

s 

8 

0.16 

0.0043 




E 

N 

w* 

c 

s 

12 

0.36 

0.0107 

LEGEND: E = essential boundary condition, N = natural boundary condition; 

L= left essential boundary condition; R= right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS = conditional stability; U = unstable 
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Table 12. SOLUTION PROCEDURES AND RESULTS USING CLASSICAL 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAINS 
TWO AND THREE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

(«;)o 


Stab. 

# of 

% 

CPU’ 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0 < x<2 

Prev. 

Midpt. 

E 

E 

(i) 

c 

s 

17 

0.38 

0.0678 

Value 

E 

N 

a) 

c 

s 

77 

0.34 

0.2262 

0 <jr<2 

Prev. 

Linear 

E 

E 

(i) 

c 

s 

17 

0.25 

0.0458 

Value 

E 

N 

(i) 

c 

s 

76 

0.29 

0.2239 

0 < jc < 2 

Avg. 

Midpt. 

E 

E 

(i) 

c 

s 

12 

0.36 

0.0423 

Value 

E 

N 

a) 

c 

s 

18 

0.34 

0.0572 

0 < jc < 2 

Avg. 

Linear 

E 

E 

a) 

c 

s 

12 

0.23 

0.0298 

Value 

E 

N 

(0 

c 

s 

18 

0.25 

0.0465 

0<jc<5 


Midpt. 

E 

E 

(i) 

c‘ 

s 

200 + 

0.60 

2.0412 


E 

N 

(i) 

c" 

s 

200 + 

0.96 

3.2696 

0 < x < 5 


Linear 

E 

E 

(i) 

c‘ 

s 

200 + 

0.43 

1.4574 


E 

N 

(i) 

c‘ 

s 

200 + 

0.65 

2.2315 

0 < *<5 

Avg. 

Midpt. 

E 

E 

(i) 

c 

s 

13 

0.60 

0.1365 

Value 

E 

N 

(i) 

c 

s 

10 

0.65 

0.0986 

0<x<5 

Avg. 

Linear 

E 

E 

0 ) 

c 

s 

13 

0.42 

0.0777 

Value 

E 

N 

(i) 

c 

s 

13 

0.42 

0.0863 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R= right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS= conditional stability; U = unstable _ 

NOTES: a - Convergence was imminent and obtained within another 50 iterations 
(1) - See equation (4.10) 
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Table 13. SOLUTION PROCEDURES AND RESULTS USING CLASSICAL 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAIN 
ONE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

W) 0 

g 

Stab. 

# of 

% 

CPU’ 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 




E 

E 

a 

c 

S 

9 

6.94 

0.7851 

0 < x< 1 

Prev. 

Midpt. 

E 

E 

M 

c 

S 

8 

6.98 

0.6507 

Value 

E 

E 

R 

c 

S 

10 

6.91 

0.8740 




E 

N 

n 

c° 

S 

17 

.... 





E 

E 

n 

c 

s 

9 

4.61 

0.4300 

0 < x < 1 

Prev. 

Linear 

E 

E 

M 

c 

s 

8 

4.66 

0.5275 

Value 

E 

E 

R 

c 

s 

10 

4.58 

0.5785 




E 

N 

a 

WBU 

s 

2004- 

.... 





E 

E 

a 

c 

s 

13 

6.91 

1.0120 

0 <x < 1 

Avg. 

Midpt. 

E 

E 

M 

c 

s 

11 

7.04 

0.9376 

Value 

E 

E 

R 

c 

s 

15 

6.90 

1.1933 




E 

N 

ES 

c° 

s 

17 

.... 





E 

E 

n 

c 

s 

12 

4.57 

0.5988 

0 < or< 1 

Avg. 

Linear 

E 

E 

M 

c 

s 

11 

4.71 

0.6376 

Value 

E 

E 

R 

c 

s 

15 

4.57 

0.8058 ! 




E 

N 

a 

C' 

s 

17 

.... 



LEGEND: E = essential boundary - condition; N = natural boundary condition; 

L= left essential boundary condition; R= right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S - unconditional stability; 
CS = conditional stability; U = unstable _ 

NOTES: a - Converges to another solution of the differential equation 

b - Convergence to another solution of the differential equation was im¬ 


minent 
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Table 14. SOLUTION PROCEDURES AND RESULTS USING CLASSICAL 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAINS 
TWO AND THREE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

K)< 


Stab. 

# of 

0 / 

VO 

CPU’ 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0 < x < 2 

Prev. 

Midpt. 

E 

E 

(l) 

C* 

S 

87 

.... 


Value 

E 

N 

(i) 

( 2 ) 

.... 

200 + 

.... 


0 < x < 2 

Prev. 

Linear 

E 

E 

(i) 

■a 

S 

95 

.... 


Value 

E 

N 

O) 

( 2 ) 

.... 

200 + 

.... 


0 < x < 2 

Avg. 

Midpt. 

E 

E 



s 

40 

.... 


Value 

E 

N 

(i) 

C‘ 

s 

42 

.... 


0 <jc< 2 

Avg. 


E 

E 

(l) 

c 

s 

35 

.... 


Value 

Linear 

E 

N 

(i) 

c 

s 

52 

—— 


0<x<5 


Midpt. 

E 

E 

(i) 

(2) 

.... 

200 + 

.... 



E 

N 

(i) 

( 2 ) 

—— 

200 + 

—— 


0<x< 5 


Linear 

E 

E 

(i) 

( 2 ) 

.... 

200 + 

.... 



E 

N 

(i) 

( 2 ) 

—. 

200 + 

—- 


0 < x < 5 

Avg. 

Midpt. 

E 

E 

(i) 

C‘ 

s 

117 

.... 


Value 

E 

N 

0) 

c' 

s 

80 

—— 


0 < x < 5 

Avg. 

Linear 

E 

E 

(i) 

c* 

s 

43 

.... 


Value 

E 

O 

(i) 

c' 

s 

105 

.... 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R * right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS = conditional stability; U = unstable 


NOTES: a - Converges to a nonsolution of the differential equation 

(1) - See equation (4.12) 

(2) - Nonconvergent solution 


Three general observations can be made from the results shown in these 
tables. The first is that, for solution procedures which yielded valid converged approxi¬ 
mations, the number of iterations to convergence was strictly a function of the initial and 
subsequent iteration strategies and independent of the interpolation strategy. Addi¬ 
tionally, the number of iterations to convergence was independent of the number of el¬ 
ements utilized. The accuracy of the solution, on the other hand, w r as strictly a function 





























































































































of the interpolation strategy and the number of degrees of freedom used in the approxi¬ 
mation. 

The second observation is that the classical linearization technique provided 
valid approximations of equation (4.1) over domain three, but was not able to solve 
equation (4.2) over domain two, despite the fact both of these equations over these re¬ 
spective domains have the same function order. The main reason for this lies in the fact 
that the solution of equation (4.2) over domain two is changing more rapidly than that 
of equation (4.1) over domain three. Thus, it appears that the classical linearization 
technique can provide valid approximations of nonlinear differential equations that have 
a function order of two, but whose rate of change over the given domain is not 'exces¬ 
sively' large. 

Last of all, a general comment on the magnitude of the average percent 
difference values in Table 13 is in order. Though these values may seem large in relation 
to the percent difference values in Table 11 and Table 12, it must be remembered that 
these are average values. In this particular situation, the values of the dependent vari¬ 
able over 0 < x < 0.2 are quite small, i.e., on the order of 0.08 and less. Thus, an ap¬ 
proximate solution at one node which is in absolute error of only 0.0023 from an exact 
solution value of 0.0100 yields a 23 percent error. Therefore, the larger average percent 
difference values stem from minor errors in the approximations at those nodes where the 
magnitude of the dependent variable is very small. The overall approximate solutions 
provided by these procedures is much better than the average percent difference values 
indicate as the actual percent difference values over a majority of the domain was less 
than 0.5 percent. The effect of each problem parameter on the performance of the var¬ 
ious solution procedures which provided convergent solutions of equations (4.1) and 
(4.2) is now examined. 

b. Boundary Conditions 

The use of an E-E boundary condition combination provided convergence 
with less iterations than the E-N combination in all instances except for the solution of 
equation (4.1) over domain three. No explanation for this behavior could be determined 
and it remains an open question requiring further investigation. The accuracy of the 
approximations within a specific combination of iteration and interpolation strategies 
was not greatly affected by the boundary conditions, as the average percent difference 
values provided by the E-E and E-N combinations were similar. 
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c. Initial Iteration Strategy 

In solving equations (4.1) and (4.2) over domain one, the initial iteration 
strategy is developed utilizing the prescribed essential boundary condition(s) as described 
in III.C.l. An initial iteration strategy based on the average value of the two essential 
boundary' conditions consistently provided convergence with slightly less iterations than 
the other two strategies. 

Over domains two and three, only convergent solutions of equation (4.1) 
were obtained; no valid approximate solutions were obtained of equation (4.2) over these 
domains. Initial iteration strategies defined by equations (4.9) and (4.10) were both 
utilized in the analysis over domain two to determine which was most effective. 
Equations (4.9) and (4.10) provided nearly the same numerical values in the approxi¬ 
mation, but the use of equation (4.10) consistently enabled the solution procedure to 
converge with less iterations. This behavior requires further research and remains an 
open question, as it was felt that equation (4.9), which accounted for the dominance of 
the u” term over the first part of the domain, should have provided better initial iteration 
values. In the solution of equation (4.2) over domains two and three, both initial iter¬ 
ation strategies led to convergent/nonconvergent nonsolutions of the differential 
equation. 

d. Subsequent Iteration Strategy 

Over domain one for both equations (4.1) and (4.2), the use of the previous 
value iteration strategy consistently provided convergence with slightly less iterations 
than the average value strategy. The one exception to this was when an E-N boundary 
condition combination was specified. In that case, both strategies yielded convergence 
with the same number of iterations. Over domains two and three, the average value 
strategy provided for convergence with significantly less iterations than the previous 
value method. This fact is especially evident in the solution of equation (4.1) over do¬ 
main three where the previous value strategy could not converge within 200 iterations 
while the average value method provided convergence in 13 iterations or less. These two 
observations tend to indicate that the classical linearization strategy provides for 
monotonic convergence when the function order is one or less and oscillatory conver¬ 
gence when the function order is greater than one. 

e. Interpolation Strategy 

The linear interpolation strategy consistently provided more accurate ap¬ 
proximations than the midpoint strategy as indicated by the lower average percent dif¬ 
ference values in Table 11 through Table 13. The main reason for this relates to the 


61 







fact that the solutions of equations (4.1) and (4.2) are polynomial functions of x. The 

r,i 

number of elements utilized in each of the solution procedures was sufficient enough to 
make the solution curve almost linear over each element. Thus, the linear approxi¬ 
mation induced less error in the evaluation of the Galerkin linearization matrix integral 
than the midpoint method. For those cases where only a few elements were utilized, the 
midpoint method occasionally provided a more accurate approximation than the linear 
strategy. But, the overall accuracy of the approximation was poor due to the decrease 
in the number of system DOF. 

/. Overall Performance 

When the function order of the differential equation solution is one or less, 
a solution procedure which utilizes a previous value iteration and linear interpolation 
strategy tends to provide the most efficient solutions based on the average percent dif¬ 
ference values in Table 11 and CPU' values in Table 13. The CPU’ values in Table 11 
should be evaluated with caution as the CPU times upon which they are based are at 
or below the accuracy level of the clock subroutine which is ± 0.03 seconds. When the 
function order is greater than one, an average value iteration and linear interpolation 
strategy combination provide the most efficient approximation as indicated by the 
CPU’ values in Table 12. 

3. Conclusions 

The classical linearization strategy can generally provide approximations of 
second order nonlinear differential equations when the function order is two or less and 
the rate of change of the dependent variable is not extremely large. Provided these 
conditions are met, an average value iteration strategy combined with a linear interpo¬ 
lation strategy should provide an efficient, valid approximation of the differential 
equation. If the function order is later found to be one or less, the use of a previous 
value iteration strategy should result in similar numerical results and converge using 
slightly less iterat ns. In either case, both the E-E and E-N boundary condition com¬ 
binations can be accommodated, although the use of an E-E combination is preferred. 

D. QUASILINEARIZATION 
1. Problem Formulation 

Quasilinearization transforms the nonlinear tf term into a linear term using the 
relation given by equation (3.4.a), where q(u ) = (u’) J and q'(u ') = 2u\ Substitution of 
these functions into equation (3.4.a) yields 
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(4.22) 


= 2 u u + (( uf — u(2u)) 

= 2 u u — (uf 

where u' is determined as outlined in III.C. Substitution of equation (4.22) into 
equations (4.1) and (4.2) yields a linear differential equation of the form 

u" ± 2u u — f(x) ± (uf xeD (4-23) 

where both ' signs are for equation (4.1), both ' + ' signs are for equation (4.2) and 
f(x) is the respective excitation function in each equation. The Galerkin FEM formu¬ 
lation process transforms equation (4.23) into 

G(G t )'u | * - f G'(G T )'rfx u ± 2 f GG*udx u = f Gf(x)dx ± f G(ufdx (4.24) 
O Jp Jp Jp Jp 

The first two terms on the left side of equation (4.24) yield B — Au where the B vector 
is present only when a natural boundary condition is prescribed. The third term on the 
left side of the equation is the linearization matrix integral which was evaluated in IV.C.l 
and yields L\ Both integrals on the right side of equation (4.24), the excitation and 
linearization vector, were evaluated in IV.B.I and yield F and F’, respectively. 

Substitution of the matrix and vector equivalents for each term into equation 
(4.24) yields a system of equations given by 

B - Au + 2L*u = F + F* (4.25.a) 

where the ± signs in equation (4.24) are incorporated into L’ and F\ Combining F and 
B to yield F m , the system of linear algebraic equations takes the final form 

(“ A + 2L' 1 )u,= F m + F*, (4.25.4) 

Equation (4.25.b) is solved iteratively for u, with both L' and F’ changing after each it¬ 
eration, where i is the iteration counter. 

2. Results 

a. General 

Forty eight different solution procedures were evaluated in solving 
equations (4.1) and (4.2) over domain one while there were twenty four procedures 
available for approximating each equation over both domains two and three. The 
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FORTRAN programs for quasilinearization are provided in Appendix G. The different 
strategies utilized in each solution procedure and the corresponding results in approxi¬ 
mating the solution of equation (4.1) are provided in Table 15 on page 65 through 
Table 19 on page 69. For equation (4.2), this information is provided in Table 20 on 
page 70 through Table 25 on page 75. The number of elements utilized in each solution 
procedure is given in Table 5 on page 47. 

Four general observations can be made about the quasilinearization strat¬ 
egy based on these results. The first is that this linearization technique provided valid 
approximations of both equations (4.1) and (4.2) over all domains, although some indi¬ 
vidual solution procedures were much more accurate than others. The difference in 
performance between the various solution procedures became noticeable as the function 
order approached three or more and the solution function gradient became large, as in 
equation (4.2) over domains two and three. 

The second point is that this linearization technique provides for conver¬ 
gence with a minimum of iterations due to its quadratic rate of convergence [Ref. 1: pp. 
38-40]. Some of the solution procedures utilized to solve equations (4.1) and (4.2) over 
domain three converged in just two iterations and thus did not yield as accurate sol¬ 
utions as was anticipated. In order to determine if these solution procedures could 
provide more accurate approximations, the convergence criterion was changed from 
.0001 to .0000001, to allow for slightly more iterations. The effect of changing the con¬ 
vergence criterion for those applicable solution procedures is noted in Table 18, 
Table 19, Table 24, and Table 25. 

The third point, which also relates to convergence, is that the number of 
iterations required for convergence is not always a function of the overall iteration 
strategy as it was for the previous two linearization strategies. When the function order 
is one or slightly over two, as in equation (4.1) over domains one and tw’o, and equation 
(4.2) over domain one, the number of iterations is dictated by the iteration strategies 
utilized. For those situations where the function order is almost three or more, the 
number of iterations required for convergence also appears to be affected by the specific 
combination of interpolation strategies used in the solution procedure. 
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Table 15. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAIN 
ONE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

("Do 


Stab. 

1 # of 

0/ 

7 0 

CPU - 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 




E 

E 

n 

C 

s 

4 

0.00 

0.0000 

0< x< 1 

Prev. 

Midpt. 

E 

E 

M 

C 

s 

4 

0.00 

0.0000 

Value 

Midpt. 

E 

E 

R 

c 

s 

4 

0.00 

0.0000 




E 

N 

n 

c 

s 

5 

0.00 

0.0000 




E 

E 

n 

c 

s 

4 

0.60 

0.0100 

0< x < 1 

Prev. 

Midpt. 

E 

E 

M 

c 

s 

4 

0.60 

0.0080 

Value 

1 ,'4-3/4 

E 

E 

R 

c 

s 

4 

0.60 

0.0100 




E 

X 

a 

c 

s 

5 

0.16 

0.0023 




E 

E 

n 

c 

s 

4 

0.32 

0.0053 

0<x< 1 

Prev. 

Midpt. 

E 

E 

M 

c 

s 

4 

0.32 

0.0064 

Value 

Linear 

E 

E 

R 

c 

s 

4 

0.32 

0.0064 




E 

E 

R 

c 

s 

5 

0.24 

0.0057 




E 

E 

m 

c 

s 

4 

0.16 

0.0021 

0 < jr < 1 

Prev. 

Linear 

E 

E 

M 

c 

s 

4 

0.16 

0.0021 

Value 

Midpt. 

E 

E 

R 

c 

s 

4 

0.16 

0.0032 




E 

X 

a 

c 

s 

5 

0.61 

0.0122 




E 

E 

a 

c 

s 

4 

0.44 

0.0059 

0 < jr < 1 

Prev. 

Linear 

E 

E 

M 

c 

s 

4 

0.44 

0.0059 

Value 

1/4-3,'4 

E 

E 

R 

c 

s 

4 

0.44 

0.0044 




E 

N 

a 

c 

s 

5 

0.54 

0.0108 




E 

E 

MM 

c 

s 

4 

0.16 

0.0027 

A x v ✓ 1 

Prev. 

Linear 

E 

E 

M 

c 

s 

4 

0.16 

0.0027 

U v X < 1 

Value 

Linear 

E 

E 

R 

c 

s 

4 

0.16 

0.0027 




E 

N 

n 

c 

s 

5 

0.02 j 

0.0004 

LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R = right essential boundary condition; 

M = average of L and R; C - convergence; D = divergence S = unconditional stability; 
CS* conditional stability; U = unstable 




































































































































































































Table 16. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAIN 
ONE (CONT.) 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. ( 

Z.'s 

K)o 

C 

1 

Stab. 

# of 

% 

CPU' 


Strat. 

Strat. 

Lt 

Rt 


Iter. 

Dif 

(sec) 




E 

E 

D 

C 

s 

5 

0.00 

0.0000 


Avg. 

Midpt. 

E 

E 

M 

c 

s 

5 


0< x< 1 

Value 

Midpt. 

E 

E 

R 

c 

s 

6 

0.00 

0.0000 




E 

N 

m 

c 

s 

5 

0.02 

0.0004 




E 

E 

n 

c 

s 

5 

0.60 

0.0120 

0< x< 1 

Avg. 

Midpt. 

E 

E 

M 

c 

s 

5 

0.60 

o.ooso 

Value 

1/4-3,'4 

E 

E 

R 

c 

s 

6 

0.61 

0.0141 




E 

N 

a 

c 

s 

5 

0.11 

0.0015 




E 

E 

a 

c 

s 

5 

0.32 

0.0053 

0<jr< 1 

Avg. 

Midpt. 

E 

E 

M 

c 

s 

5 

0.32 

0.0064 

Value 

Linear 

E 

E 

R 

c 

s 

6 

0.32 

0.0075 




E 

N 

a 

c 

s 

5 

0.27 

0.0053 




E 

E 

n 

c 

s 

5 

0.16 

0.0027 

0<jr< 1 

Avg. 

Linear 

E 

E 

M 

c 

s 

5 

0.16 

0.0032 

Value 

Midpt. 

E 

E 

R 

c 

s 

6 

0.16 

0.0038 




E 

wa 

a 

c 

s 

5 

0.60 

0.0099 




E 

E 

n 

c 

s 

5 

0.45 

0.0089 

0< x< 1 

Avg. 

Linear 

E 

E 

M 

c 

s 

5 

0.45 

0.0089 

Value 

1/4-3/4 

E 

E 

R 

c 

s 

6 

0.44 

0.0118 




E 

N 

n 

c 

s 

5 

0.53 

0.0070 




E 

E 

a 

c 

s 

5 

0.16 

0.0037 

0<x< 1 

Avg. 

Linear 

E 

E 

M 

c 

s 

5 

0.16 

0.0038 

Value 

Linear 

E 

E 

R 

c 

s 

6 

0.16 

0.0037 




E 

N 

a 

c 

s 

5 

0.36 

0.0084 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R= right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS » conditional stability; U = unstable_ 
























































































































































































































Table 17. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAIN 
TWO 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

K)o 


Stab. 

#of 

% 

CPU’ 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0 < or < 2 

Prev. 

Midpt. 

E 

E 

(1) 

C 

s 

4 

0.00 

0.0001 

Value 

Midpt. 

E 

N 

(1) 

C 

s 

4 

0.00 

0.0001 

0< jc< 2 

Prev. 

Midpt. 

E 

E 

(1) 

c 

s 

4 

0.97 

0.0547 

Value 

1/4-3/4 

E 

N 

(1) 

c 

s 

4 

0.89 

0.0475 

0< x< 2 

Prev. 

Midpt. 

E 

E 

(1) 

c 

s 

4 

0.52 

0.0291 

Value 

Linear 

E 

N 

(1) 

c 

s 

5 

0.47 

0.0249 

0< x < 2 

Prev. 

Linear 

E 

E 

(1) 

c 

s 

4 

0.26 

0.0130 

Value 

Midpt. 

E 

Q 

(1) 

c 

s 

4 

0.22 

0.0119 

0<x< 2 

Prev. 

Linear 

E 

E 

(1) 

c 

s 

4 

0.71 

0.0424 

Value 

1 ,'4-3/4 

E 

O 

(1) 

c 

s 

4 

0.69 

0.0411 

0 < x < 2 

Prev. 

Linear 

E 

E 

(1) 

c 

s 

4 

0.26 

0.0137 

Value 

Linear 

E 

N 

(1) 

c 

s 

4 

0.27 

0.0153 

0 < or < 2 

Avg. 

Midpt. 

E 

E 

(1) 

c 

s 

6 

0.00 

0.0001 

Value 

Midpt. 

E 

N 

(1) 

c 

s 

5 

0.00 

0.0002 

0<jr<2 

Avg. 

Midpt. 

E 

E 

(1) 

c 

s 

5 

0.97 

0.0610 

Value 

1/4-3,'4 

E 

N 

(1) 

c 

s 

5 

0.89 

0.0594 

0 < x < 2 

Avg. 

Midpt. 

E 

E 

(1) 

c 

s 

5 

0.51 

0.0308 

Value 

Linear 

E 

N 

(1) 

c 

s 

5 

0.47 

0.0326 

0 < or < 2 

Avg. 

Linear 

E 

E 

(1) 

c 

s 

6 

0.26 

0.0198 

Value 

Midpt. 

E 

N 

(1) 

c 

s 

6 

0.22 

0.0179 

0<x<2 

Avg. 

Linear 

E 

E 

(1) 

c 

s 

5 

0.71 

0.0448 

Value 

1 ,'4-3/4 

E 

N 

(1) 

c 

s 

5 

0.69 

0.0503 

0< jc< 2 

Avg. 

Linear 

E 

E 

(1) 

c 

s 

5 

0.26 

0.0163 

Value 

Linear 

E 

N 

(1) 

c 

s 

5 

0.27 

0.0180 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R = right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS conditional stability; U - unstable _ 

NOTES: (1) - See equation (4.10)_ 
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Table 18. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAIN 
THREE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

K) 0 


Stab. 

# of 

% 

CPU’ 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0<jc<5 

Prev. 

Midpt. 

E 

E 

(i) 

C° 

S 

5 

0.00 

0.0002 

Value 

Midpt. 

E 

N 

(i) 

C‘ 

S 

6 

0.04 

0.0046 

0<jc<5 

Prev. 

Midpt. 

E 

E 

U) 

c 

S 

3 

1.43 

0.0948 

Value 

1;4-3/4 

E 

N 

(i) 

c 

S 

7 

1.55 

0.2216 

0 < x < 5 

Prev. 

Midpt. 

E 

E 

(i) 

■a 

S 

2 

0.53 

0.0230 

Value 

Linear 

E 

N 

O) 

c 

S 

5 

0.83 

0.0911 

0 < x < 5 

Prev. 

Linear 

E 

E 

(i) 

c" 

s 

5 

0.38 

0.0352 

Value 

Midpt. 

E 

N 

(l) 

c 

s 

5 

0.43 

0.040S 

0 <x < 5 

Prev. 

Linear 

E 

E 

(i) 

c* 

s 

2 

0.80 

0.0400 

Value 

1'4-3/4 

E 

N 

(i) 

c 

s 

3 

1.08 

0.0717 

0< x< 5 

Prev. 

mm 

E 

E 

(i) 

C° 

s 

2 

0.15 

0.0063 

Value 

■sal 

E 

N 

(i) 

El 

s 

2 

0.15 

0.0068 


LEGEND: E * essential boundary condition; N = natural boundary condition; 

L= left essential boundary condition; R = right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S= unconditional stability; 
CS = conditional stability; U = unstable _ 

NOTES: a - Convergence criterion changed to 0.0000001 in order to allow for ad¬ 
ditional iterations and a more efficient approximation 

b - Convergence criterion kept at .0001 as use of tighter criterion in note 
a leads to a less efficient or nonconvergent approximation. 

(1) - See equation (4.10) 
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Table 19. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.1) OVER DOMAIN 
THREE (CONT.) 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

("Do 


Stab. 

# of 

% 

CPU* 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0 < x < 5 

Avg. 

Midpt. 

E 

E 

(i) 

C° 

S 

5 

0.00 

0.0001 

Value 

Midpt. 

E 

N 

(i) 

C‘ 

s 

9 

0.04 

0.0069 

0< x< 5 

Avg. 

Midpt. 

E 

E 

(i) 

W3M 

s 

2 

2.95 

0.1276 

Value 

1/4-3/4 

E 

N 

(i) 

C 

s 

9 

1.55 

0.2782 

0 < x < 5 

Avg. 

Midpt. 

E 

E 

(i) 

c* 

s 

7 

0.76 

0.1040 

Value 

Linear 

E 

N 

(i) 

c 

s 

6 

0.83 

0.1049 

0 < x < 5 

Avg. 

Linear 

E 

E 

(i) 

c° 

s 

7 

0.38 

0.0528 

Value 

Midpt. 

E 

N 

(i) 

c 

s 

5 

0.42 

0.0422 

0<x< 5 

Avg. 

Linear 

E 

E 

(i) 

c‘ 

s 

8 

1.05 

0.1644 

Value 

1 ,'4-3/4 

E 

N 

(i) 

c 

s 

4 

1.07 

0.0818 

0 < x < 5 

Avg. 

'-X 

E 

E 

(i) 

c 

s 

2 

3.35 

0.1559 

Value 


E 

N 

(i) 

c‘ 

s 

6 

0.42 

0.0519 


LEGEND: E = essential boundary condition; N = natural boundary condition; 
L=left essential boundary condition; R= right essential boundary condition; 

M = average of L and R; C = convergence; D » divergence S = unconditional stability; 
CS = conditional stability; U = unstable _ 

NOTES: a - Convergence criterion changed to 0.0000001 in order to allow for ad¬ 
ditional iterations and a more efficient approximation 

b • Convergence criterion kept at .0001 as use of tighter criterion in note 
a leads to a less efficient or nonconvergent approximation. 

c - Tightening of convergence criterion had no effect 
_(1) - See equation (4.10)_ 
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Table 20. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAIN 
ONE 



Solution Procedure 

Results 

Iter. 

Strat. 

Interp. 

Strat. 

B. C.'s 

K)o 


Stab. 

# of 
Iter. 

% 

Dif 

CPU* 

(sec) 

Lt 

Rt 

0< x< 1 

Prev. 

Value 

Midpt. 

Midpt. 

E 

Q 

a 

C 

s 

5 

1.05 

0.0663 

E 

E 

M 

c° 

s 

11 

.... 


E 

E 

R 

c 

s 

5 

.... 


E 

N 

m 

c° 

s 

6 

.... 


0<jr< I 

Prev. 

Value 

Midpt. 
1/4-3,'4 

E 

E 

n 

c 

s 

5 

16.1 

1.0700 

E 

E 

M 

C‘ 

s 

12 

.... 


E 

E 

R 

- ,0 

s 

5 

.... 


E 

N 

m 

c' 

s 

6 

.... 


0<x< 1 

Prev. 

Value 

Midpt. 

Linear 

E 

E 

b 

c 

s 

5 

9.28 

0.6487 

E 

E 

M 

c‘ 

s 

12 

.... 


E 

E 

R 

c° 

s 

5 

.... 


E 

N 

B 

c" 

s 

6 

.... 


0 < jr < 1 

Prev. 

Value 

Linear 

Midpt. 

E 

E 

b 

c 

s 

5 

3.61 

0.2163 

E 

E 

M 

c° 

s 

11 

.... 


E 

E 

R 

C' 

s 

5 

.... 


E 

N 

B 

C‘ 

s 

6 

.... 


0 <x< 1 

Prev. 

Value 

Linear 
1/4-3/4 

E 

E 

B 

c 

s 

5 

11.4 

0.6809 

E 

E 

M 

c* 

s 

12 

.... 


E 

E 

R 

c' 

s 

5 

.... 


E 

N 

B 

c' 

s 

6 

.... 


0 < jc < 1 

Prev. 

Value 

Linear 

Linear 

E 

E 

B 

c 

s 

5 

4.64 

0.3087 

E 

E 

M 

c” 

s 

12 

.... 


E 

E 

R 

c* 

s 

5 

.... 

[cm 

E 

N 

B 

c* 

s 

6 

.... 


LEGEND: E- essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R = right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS = conditional stability; U = unstable 

NOTES: a - Converges to another solution of the differential equation 
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Table 21. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAIN 
_ ONE (CONT.) _ 


Domain 

Solution Procedure 

Results 

Iter. 

Strat. 

Interp. 

Strat. 

B. C.'s 

(«;)„ 


Stab. 

# Of 

Iter. 

% 

Dif 

CPU- 

(sec) 

Lt 

Rt 

0 < jc < 1 

Avg. 

Value 

Midpt. 

Midpt. 

E 

E 

B 

c 

s 

6 

1.07 

0.0852 

E 

E 

M 

c° 

s 

15 

.... 


E 

E 

R 

c' 

s 

7 

.... 


E 

N 

B 

c° 

s 

7 

.... 


0< x< I 

Avg. 

Value 

Midpt. 

I/4-3/4 

E 

E 

B 

c 

s 

6 

16.1 

1.3399 

E 

E 

M 

c° 

s 

16 

.... 


E 

E 

R 

c° 

s 

7 

.... 


E 

N 

B 

c° 

s 

7 

.... 


0< x< 1 

Avg. 

Value 

Midpt. 

Linear 

E 

E 

B 

c 

s 

6 

9.30 

0.7422 

E 

E 

M 

c° 

s 

15 

.... 


E 

E 

R 

C° 

s 

7 

.... 


E 

N 

B 

c 

s 

7 

.... 


0 < jc < 1 

Avg. 

Value 

Linear 

Midpt. 

E 

E 

B 

c 

s 

7 

3.61 

0.3486 

E 

E 

M 

C‘ 

s 

15 

.... 


E 

E 

R 

C 

s 

6 

.... 


E 

N 

B 

c" 

s 

7 

.... 


0 <jC < 1 

Avg. 

Value 

Linear 

1,'4-3/4 

E 

E 

B 

c 

s 

7 

11.4 

0.9091 

E 

E 

M 

c° 

s 

16 

.... 


B 

E 

R 

C' 

s 

7 

.... 


E 

X 

B 

■a 

s 

7 

.... 


0 < x < 1 

Avg. 

Value 

■ 

E 

E 

B 

c 

s 

7 

4.60 

0.4440 

E 

E 

M 

c* 

s 

15 

.... 


E 

E 

R 

c 

s 

7 

.... 


E 

X 

B 

c" 

s 

7 

.... 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R = right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS = conditional stability; U = unstable 

NOTES: a - Converges to another solution of the differential equation 
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Table 22. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAIN 
TWO 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. < 

Z.'s 

K) 0 


Stab. 

# Of 

% 

CPU - 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0 < x < 2 

Prev. 

Midpt. 

E 

E 

(i) 

c 

S 

4 

4.62' 

1.1534 

Value 

Midpt. 

E 

N 

(l) 

C 

S 

3 

0.85 

0.1777 

0<x<2 

Prev. 

Midpt. 

E 

E 

(i) 

c 


KM 

56.6' 

15.252 

Value 

1 .'4-3/4 

E 

N 

U) 

—— 

BM 

200 + 

d 


0 < x < 2 

Prev. 

Midpt. 

E 

E 

0) 

c 

s 

4 

32.0' 

9.0524 

Value 

Linear 

E 

\ 

(1) 

.... 

mssm 

200 + 

d 


0 < x < 2 

Prev. 

Linear 

E 

E 

(1) 

c 

s 

3 

9.68 / 

2.0618 

Value 

Midpt. 

E 

N 

(1) 

c 

s 

5 

EM 

92.669 

0 < x < 2 

Prev. 

Linear 

E 

E 

(1) 

c 

ESI 

4 

38.9' 

10.626 

Value 

1/4-3,'4 

E 

N 

(1) 

c 

s 

5 

188' | 

64.094 

0<x< 2 

Prev. 

Linear 

E 

E 

(1) 

c 

s 

4 

Tiftii 

4.3182 

Value 

Linear 

E 

N 

(1) 

c 

s 

3 

04 

OO 

o 

0.6880 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L=left essential boundary condition; R 3 * right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 

CS = conditional stability; U = unstable _ 

NOTES: a - Majority of error occurs over 0 < x < 0.2 

b - Divergence results when certain number of elements are utilized 
c - Majority of error occurs over 0 < x < 0.6 
d - Provides a reasonable approximation over 1 < x < 2 
e - Majority of error occurs over 0 < x < 0.5 
f - Majority of error occurs over 0 < x < 0.3 


(1) - See equation (4.11) 





















































































































Table 23. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAIN 
TWO (CONT.) _ 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

K)o 


— 

Stab. 

# of 

% 

CPU’ 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0 < x < 2 

Avg. 

Midpt. 

E 

E 

(i) 

C 

s 

4 

4.70° 

1.2206 

Value 

Midpt. 

E 

N 

0) 

C 

s 

4 

0.92 

0.2580 

0<x<2 

Avg. 

Midpt. 

E 

E 

(i) 

c 

s 

4 

Hindi 

14.952 

Value 

1/4-34 

E 

N 

(i) 

c 

cs c 

129 

EES 

4089.7 

0<x<2 

Avg. 

Midpt. 

E 

E 

(i) 

c 

s 

4 

31.5' 

9.1124 

Value 

Linear 

E 

N 

(i) 

c 

s 

21 

KEEH 

811.74 

0 <x< 2 

Avg. 

Linear 

E 

E 

(i) 

c 

s 

4 


2.4966 

Value 

Midpt. 

E 

N 

(i) 

c 

s 

5 

Esa 

88.739 

0 <x <2 

Avg. 

Linear 

E 

E 

(i) 

c 

s 

4 

38.8' 

11.233 

Value 

1/4-3/4 

E 

N 

(i) 

c 

s 

5 

TEH 

68.267 

0 <x <2 

Avg. 


E 

E 

(i) 

c 

s 

4 

16.0' 

4.0571 

Value 


E 

N 

0) 

c 

s 

4 

3.12° 

0.8615 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R= right essential boundary condition; 

M * average of L and R; C = convergence; D = divergence S = unconditional stability; 

CS = conditional stability; U = unstable _ 

NOTES: a - Majority of error occurs over 0 < x < 0.2 
b - Majority of error occurs over 0 < x < 0.6 
c - Divergence results when certain number of elements are utilized 
d - Provides a reasonable approximation over 1 < x < 2 
e - Majority of error occurs over 0 < x < 0.5 
/- Majority of error occurs over 0 < x < 0.3 

_(1) - See equation (4.11)_ 
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Table 24. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAIN 
THREE 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C. s 

«)o 


Stab. 

# Of 

% 

CPU’ 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0 < x < 5 

Prev. 

Midpt. 

E 

E 

(1) 

c° 

S 

5 

0.16 

0.0875 

Value 

Midpt. 

E 

N 

(1) 

c‘ 

S 

11 

0.38 

0.4189 

0 < jc < 5 

Prev. 

Midpt. 

E 

E 

(1) 

D 

U 


.... 


Value 

1/4-3/4 

E 

N 

(1) 

D 

um 


.... 


0 < x < 5 

Prev. 

Midpt. 

E 

E 

(1) 

D 

u 




Value 

Linear 

E 

N 

(1) 

D 

u 


.... 


0 < x < 5 

Prev. 

Linear 

E 

E 

(I) 

El 

CS c 

2 


636.11 

Value 

Midpt. 

E 

N 

(1) 

D 

mm 


.... 


0< x< 5 

Prev. 

Linear 

E 

E 

(1) 

D 

u 


.... 


Value 

1/4-3,'4 

E 

N 

(1) 

D 

um 


.... 


0 < jr < 5 

Prev. 

Linear 

E 

E 

mi 


s 

5 

0.88 

0.5378 

Value 

Linear 

E 

N 

m 

S3 

s 

3 

7.48' 

2.5872 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R = right essential boundary condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS = conditional stability; L' = unstable _ 

NOTES: a - Convergence criterion changed to 0.0000001 in order to allow for ad¬ 
ditional iterations and a more efficient approximation 

b - Convergence criterion kept at .0001 as use of tighter criterion in note 
a leads to a less efficient or nonconvergent approximation. 

:.c - Divergence results when certain number of elements are used 
d - Provided an adequate approximation over 2.5 < x < 5.0 
e - Majority of error occurs over 0 < x < 0.5 where the values of u < 2.0 
_(1) - See equation (4.12) 
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Table 25. SOLUTION PROCEDURES AND RESULTS USING QUASI¬ 
LINEARIZATION TO SOLVE EQUATION (4.2) OVER DOMAIN 
THREE (CONT.) 



Solution Procedure 

Results 

Domain 

Iter. 

Interp. 

B. C.'s 

«)o 


Stab. 

# of 

% 

CPU• 


Strat. 

Strat. 

Lt 

Rt 

Iter. 

Dif 

(sec) 

0< x< 5 

Avg. 

Midpt. 

E 

E 

(1) 

C° 

S 

4 

11.1* 

5.2562 

Value 

Midpt. 

E 

N 

(1) 

c° 

S 

4 

grra 

4.8829 

0 < jc < 5 

Avg. 

Midpt. 

E 

E 

(1) 

D 

U 


—— 


Value 

1,'4-3/4 

E 

N 

(1) 

D 

U 


.... 


0< x< 5 

Avg. 

Midpt. 

E 

E 

mi 

D 

mm 


.... 


Value 

Linear 

E 

O 

(i) 

D 

■n 


.... 


0<x<5 

Avg. 

Linear 

E 


(i) 

C c 

s 

2 

DEED 

262.37 

Value 

Midpt. 

E 

N 

(i) 

D 

u 




0<x<5 

Avg. 

Linear 

E 

E 

m 

c c 

s 

2 

672' 

149.82 

Value 

1/4-3/4 

E 

N 

(l) 

D 

u 


— — 


0<x<5 

Avg. 


E 

E 

(i) 

C' 

s 

7 

1.69 

1.4075 

Value 


E 

N 

(i) 

C‘ 

s 

31 

1.46 

5.2343 


LEGEND: E = essential boundary condition; N = natural boundary condition; 

L = left essential boundary condition; R = right essential boundary- condition; 

M = average of L and R; C = convergence; D = divergence S = unconditional stability; 
CS = conditional stability; U = unstable _ 

NOTES: a - Convergence criterion changed to 0.0000001 in order to allow for ad¬ 
ditional iterations and a more efficient approximation 

b - Majority of error occurs over 0 < x < 0.5 where the values of u < 2.0 
c - Convergence criterion kept at .0001 as use of tighter criterion in note 
a leads to a less efficient or nonconvergent approximation. 

d - Provided an adequate approximation over 2.0 < x < 5.0 
e - Provided an adequate approximation over 1.5 < x < 5.0 
(1) - See equation (4.12) 


Last of all, a general comment is required with respect to the average per¬ 
cent difference values in Table 22 through Table 25. As previously noted in IV.C.3, the 
value of the dependent variable for equation (4.2) is very small over the first part of the 
domain. Thus, errors in the approximation which are on the same order of magnitude 
as the value of the dependent variable result in large percent difference values. There¬ 
fore, average percent difference values of three or more are amplified with a superscript 




































































































which advises the reader as to the true accuracy of the approximation. The term 'a 
majority of error' in each of these accompanying notes means that the percent cfltference 
at each node is 10 percent or more over the indicated domain. The effect of each pa¬ 
rameter on the overall performance of the various solution procedures is now evaluated. 

b. Boundary Conditions 

Except in a few isolated instances, the accuracy and number of iterations to 
convergence for a specific combination of iteration and interpolation strategies w r as un¬ 
affected by the boundary condition combination utilized. The one exception to this is 
shown in Table 20 and Table 21 where the use of a E-N combination caused the sol¬ 
ution procedure to converge to a second solution of the differential equation. But, when 
the initial iteration strategy was changed from using the value of the left boundary con¬ 
dition to the strategy defined by equation (4.12), valid approximations of the solution 
u = 10jc 3 were obtained by all procedures using an E-N combination. Thus, it appears 
that this linearization technique is quite favorable to both E-E and E-N type boundary 
value problems, provided that a valid initial iteration strategy is utilized. 

c. Initial Iteration Strategy 

Over domain one for equation (4.1), use of the three initial iteration strate¬ 
gies described in III.C.l provided convergence with nearly the same number of iter¬ 
ations. Over the same domain for equation (4.2), only the use cf the left essential 
boundary condition as values for the initial iteration strategy yielded a convergent ap¬ 
proximation of the original exact solution. The other two strategies, when utilized with 
an E-E boundary condition combination, converged to a second solution of equation 
(4.2). The reason that the use of the left essential boundary condition for the initial it¬ 
eration strategy outperformed the other two strategies is because it has a value of zero. 
With («;)„ = 0, the first iteration solves the differential equation neglecting the effect of 
the nonlinear term, as F* and L’ are zero. Thus, the values of («’), utilized in the next 
iteration are very close to values of (u’) 0 that would have been obtained by neglecting the 
u 1 term and integrating the differential equation twice with the imposition of boundary 
conditions. From Figure 17 on page 77, it can be seen that the u" term dominates over 
a majority of the domain. Thus, if equation (4.2) had been analyzed as order two over 
this domain and the initial iteration strategy developed by neglecting the u 2 term, similar 
results using one less iteration would probably have been obtained. 

For equation (4.1) over domain two, the initial iteration strategies defined 
by equations (4.9) and (4.10) were both used in the analysis to determine which was most 
effective. As in the classical linearization results, they both provided almost identical 
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Figure 17. Dominance of Terms in Equation (4.2) Over Domain One 
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numerical approximations, and the use of equation (4.10) again resulted in convergence 
with less iterations. Likewise, equations (4.11) and (4.12) were both utilized as initial 
iteration strategies for equation (4.2) over domain two. They both provided the same 
accuracy in the approximation, but the use of equation (4.11), which divides the 
dominance equally over the domain, converged using less iterations as was originally 
expected. The use of equations (4.10) and (4.12) for the initial iteration strategies of 
equations (4.1) and (4.2), respectively, over domain three, lead to accurate approxi¬ 
mations for those solution procedures which did converge. 

d. Subsequent Iteration Strategy 

The use of either the previous value or average value iteration strategy, in 
general, had no effect on the number of iterations required for convergence. This is most 
likely due to the quadratic rate of convergence guaranteed by the quasilinearization 
method, as previously mentioned in IV.D.2.a. 

e. Interpolation Strategy 

The quasilinearization method requires the use of two interpolation strate¬ 
gies; one for the linearization vector and one for the linearization matrix. In the in¬ 
terpolation strategy column of Table 15 through Table 25, the upper strategy is for the 
linearization matrix and the lower one is for the linearization vector. A general trend in 
th iccuracy provided by the various combinations of interpolation strategies is evident. 
The different combinations are ranked from least to most accurate in the following list, 
where the first strategy indicated is for the linearization matrix and the second is for the 
linearization vector. 

• Midpoint; 1/4-3/4 

• Linear; 1/4-3/4 

• Midpoint; Linear 

• Linear;Midpoint 

• Linear;Linear 

• Midpoint; Midpoint 

Two conclusions can be drawn from the above list. First is that the most 
accurate interpolation strategies utilize the same interpolation technique for both 
linearization integrals. Thus, if a 1/4-3,'4 interpolation strategy for the linearization 
matrix had been developed, there is a good possibility that an overall l/4-3/4;l/4-3/4 
interpolation strategy would have provided accurate approximations. Second, the least 




refined interpolation strategy, namely midpoint;midpoint, provides the most accurate 
approximations. This result is not very surprising, as in many situations, the simplest 
method provides the best results. 

f Overall Performance 

In almost all cases, the solution procedure consisting of a previous iteration 
and a midpoint;midpoint interpolation strategy provided the most efficient approxi¬ 
mations regardless of the function order of the equation or the boundary conditions 
imposed. The only case where this procedure faltered slightly was in approximating 
equation (4.2) over domain two. But, as shown in note a of Table 22 and Table 23, it 
only had a problem approximating the solution over that part of the domain where 
ucO.l. The solution procedure utilizing a previous value iteration and a linearjlinear 
interpolation strategy u*as not quite as efficient as the above solution procedure, but 
performed in an acceptable manner. 

3. Conclusions 

Quasilinearization provides a viable method of approximating nonlinear differ¬ 
ential equations that contain the u 2 term, regardless of the function order of the equation 
and the nature of the boundary conditions imposed. The use of a previous value iter¬ 
ation and either a midpoint;midpoint or linear;linear interpolation strategy should pro¬ 
vide an accurate approximation with a minimum number of iterations, provided that the 
initial iteration strategy is adeptly chosen. The actual shape of the solution curve and 
the discretization invoked, i.e., the number of elements, are the two factors most likely 
to determine which interpolation technique provides the more accurate approximation. 

E. FINAL REMARKS 

An overall solution procedure involving quasilinearization combined with a previous 
value iteration and either a midpoint;midpoint or linearjlinear interpolation strategy 
provides excellent approximations of second order, nonlinear, one dimensional, differ¬ 
ential equations which contain the u 2 nonlinear term. As the u 2 term has a more non¬ 
linear nature than some of the other nonlinear terms encountered, i.e., u"u, (u') : , etc.; it 
is felt that this solution procedure should provide viable approximations of many non¬ 
linear, second order differential equations. It cannot be overemphasized that the success 
or failure of the above solution procedure depends greatly on the initial iteration strategy 
developed. Thus, utilization of this solution procedure requires that the user have an 
in-depth understanding of the physics involved in the system being analyzed. 
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This research has provided a fundamental baseline for the future investigation of 
techniques for solving nonlinear differential equations. The following steps provide a 
logical progression for determining the actual capabilities of a Galerkin FEM solution 
procedure utilizing quasilinearization. 

• Conduct an analysis of second order, one dimensional, nonlinear, differential 
equations which contain nonlinear terms other than u 2 . These nonlinear equations 
should be of an engineering nature for which experimental data exists to allow for 
a confirmation of the results developed by the mathematical model. 

• Investigate the ability of this solution method to solve one dimensional, nonlinear, 
fourth order differential equations. This requires some modification of the in¬ 
terpolation strategies as the Galerkin FEM must utilize cubic shape functions for 
developing a fourth order differential equation approximation. 

• Extrapolate the concepts and principals developed by this and future research to 
the solution of two dimensional, second and fourth order, nonlinear differential 
equations. 
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APPENDIX A. FORCING FUNCTION FORMULATION STRATEGIES 

On an elemental level, f*G(6x)dx . becomes /''g(6(a, + £))d£ where 

• £ is the local element coordinate, 0 < f <, I, 

• a, is the sum of all element lengths prior to the element being evaluated. For equal 
length elements, a, = (/ — 1)/, where i is the element number. 

Midpoint Lumped Approximation 

This method evaluates 6(a + £) at the midpoint of the element, £ = -y and brings it 
outside the integral as a constant yielding 





£ 



4 


= 34 


a + - 2‘ 
a +~2 


Quarter/Three Quarter Lumped Approximation 

This method takes 6(a + £) inside the shape function vector yielding 


(- 4 . 1 ) 


n. 


f = 


(6(« + {)/l-|- 

<6( “ + « >(f) 

4 




{A.2.a) 


34 


The first 6(a + £) term is evaluated at £ = - j and the second term at f = -r 1 -, yielding 
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(A.2.b) 


Consistent Evaluation 

This method calculates the exact value of the Galerkin excitation integral yielding 


f * 



(« + m 



.-tiL+t-il 

/ s / 

‘e ‘e 

«L + i1 





21, 2 3 l e 

«<L jl 

21, * V, o 



3a/ e + / e 2 
3«/ e + 2/j 


M-3) 
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APPENDIX B. PROGRAM LISTINGS AND RESULTS FOR THE 

LINEAR APPLICATION OF THE GALERKIN FEM 

c 
c 
c 
c 
c 
c 

110 COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),C00RD(100), 

:ELEN,IC0RR(100,2),NEL,NSNP 

C INPUT NUMBER OF ELEMENTS AND TOTAL LENGTH OF DOMAIN 

120 PRINT*,'INPUT NUMBER OF ELEMENTS DESIRED AND TOTAL LENGTH. ' 

130 READ(6,*) NEL.TLEN 

C CALCULATE NUMBER OF NODAL POINTS AND DEFINE LEFT BOUNDARY OF 
C DOMAIN 

135 NSNP=NEL+1 
150 COORDC1)=0. 

C DETERMINE ELEMENT LENGTH, LOCAL TO GLOBAL NODAL POINT 
C CORRESPONDENCE, AND X-COORDINATE OF EACH NODAL POINT 

155 ELEN=TLEN/FLOAT(NEL) 

160 DO 169 IEL=1,NEL 

162 ICORR(IEL,1)=IEL 

163 ICORR(IEL,2)=IEL+1 

164 COORD(IEL+l)=COORD(IEL)+ELEN 

169 CONTINUE 

C CALL SUBROUTINE SYM1A TO DETERMINE A MATRIX AND F VECTOR AND SOLVE 
C THE LINEAR SYSTEM OF EQUATIONS AU = F 

170 CALL SYM1A 

C CALL SUBROUTINE UX2EXT TO DETERMINE EXACT SOLUTION 

190 CALL UX2EXT 

C CALL SUBROUTINE OUTLIN TO OUTPUT RESULTS 

200 CALL OUTLIN 

210 END 


PROGRAM LIN1 


* THIS PROGRAM SOLVES THE DIFFERENTIAL EQUATION U" 

* U(0)=0; U'(2)=4 WITH UEXACT=X**2. 


= 2., 0<X<2 
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I 


c 

c 

c 

c 

c 

c 


* SUBROUTINE SYM1A >, t 

it 

* THIS SUBROUTINE COMPUTES THE A MATRIX AND F VECTOR FOR MAIN 

* PROGRAM LIN1 AND SOLVES THE LINEAR SET OF EQUATIONS AU = F. 


* 

* 

* 

* 

v* 


100 SUBROUTINE SYM1A 


110 COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),COORD(100), 
: ELEN,ICORR(100,2),NEL,NSNP 

120 DIMENSION AE(2,2), FS1E(2), FS2E(2), VKAREA(40600) 

C ZERO OUT A MATRIX AND F VECTOR 

140 DO 210 IZ = l.NSNP 

150 FT(IZ) = 0. 

160 DO 200 JZ = l.NSNP 

170 A(IZ,JZ) = 0. 

200 CONTINUE 

210 CONTINUE 
213 ALPHA=0. 

C ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 
215 DO 375 IEL=1,NEL 

C DETERMINE ELEMENTAL A MATRIX AND ELEMENTAL F VECTOR 

220 AE(1,1)=1./ELEN 

230 AE(1,2)=(-1./ELEN) 

240 AE(2,1)=AE(1,2) 

250 AE(2,2)=AE(1,1) 

260 FS1E(1)=ELEN 

265 FS1E(2)=ELEN 

C DISTRIBUTE AE MATRICES AND FE VECTORS INTO SYSTEM A MATRIX AND F 
C VECTOR 

300 DO 370 11=1,2 

310 DO 350 JJ=1,2 

320 IN=ICORR(IEL,II) 

330 JN=ICORR(IEL,JJ) 

340 A(IN,JN)=A(IN,JN) - AE(II,JJ) 

350 CONTINUE 

360 FT(IN)=FS1E(II) + FT(IN) 

370 CONTINUE 

372 ALPHA=ALPHA + ELEN 

375 CONTINUE 

C IMPOSE KINEMATIC AND NATURAL BOUNDARY CONDITIONS 

376 A(1,1)=1. 

377 A(1,2)=0. 

378 FT(1)=0. 
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379 FT(NSNP)=FT(NSNP)-4. 

380 M=1 

390 IDGT=3 
400 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE AU = F 

410 CALL LEQT2F(A,M,NSNP,IQ,FT,IDGT,WKAREA,IER) 
420 DO 440 NEW=1,NSNP 
430 U(NEW)=FT(NEW) 

440 CONTINUE 
450 RETURN 
460 END 
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* THIS SUBROUTINE COMPUTES THE VALUE OF U=X**2 AT THE SPECIFIED * 

* NODAL POINTS FOR MAIN PROGRAM LIN1. _ _ * ****** ** 


100 SUBROUTINE UX2EXT 

110 COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),COORD(100), 
: ELEN,ICORRC100,2),NEL,NSNP 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) = COORD(NN)**2 

150 CONTINUE 
160 RETURN 
170 END 
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SOLUTION OF U" = 2. USING CONSISTENT FORCING FUNCTION 


X-COORD 

U EXACT 

U FEM 

% DIFF 

0. 000 

0.0000 

0.0000 

0.0 

0. 200 

0.0400 

0.0400 

0. 0 

0. 400 

0.1600 

0.1600 

0. 0 

0. 600 

0.3600 

0.3600 

0. 0 

0. 800 

0.6400 

0. 6400 

0. 0 

1. 000 

1. 0000 

1. 0000 

0. 0 

1. 200 

1. 4400 

1. 4400 

0. 0 

1.400 

1. 9600 

1. 9600 

0. 0 

1. 600 

2.5600 

2.5600 

0. 0 

1. 800 

3. 2400 

3. 2400 

0. 0 

2. 000 

4.0000 

4.0000 

0. 0 
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PROGRAM LIN2 


C 

C 

C 

C 

C 

C 


* 

★ 


* THIS PROGRAM SOLVES THE DIFFERENTIAL EQUATION U" = 6X, 0<X<2 

* U(0)=0; U'(2)=12 WITH UEXACT=U**3. 


* 

* 

* 

* 


110 COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),COORD(100), 
: ELEN,ICORR(100,2),NEL,NSNP 

C INPUT NUMBER OF ELEMENTS AND TOTAL LENGTH OF DOMAIN 

120 PRINT*,'INPUT NUMBER OF ELEMENTS DESIRED AND TOTAL LENGTH. ' 

130 READ(6,*) NEL,TLEN 

C CALCULATE NUMBER OF NODAL POINTS AND DEFINE LEFT BOUNDARY OF 
C DOMAIN 


135 NSNP=NEL+1 
150 C00RD(1)=0. 

C DETERMINE ELEMENT LENGTH, LOCAL TO GLOBAL NODAL POINT 
C CORRESPONDENCE, AND X-COORDINATE OF EACH NODAL POINT 

155 ELEN=TLEN/FLOAT(NEL) 

160 DO 169 IEL=1,NEL 

162 ICORR(IEL,1)=IEL 

163 IC0RR(IEL,2)=IEL+1 

164 COORD(IEL+1)=COORD(IEL)+ELEN 

169 CONTINUE 

C CALL SUBROUTINE SYM2A TO DETERMINE A MATRIX AND F VECTOR AND SOLVE 
C THE LINEAR SYSTEM OF EQUATIONS AU = F 

170 CALL SYM2A 


C CALL SUBROUTINE UX3EXT TO DETERMINE EXACT SOLUTION 


190 CALL UX3EXT 

C CALL SUBROUTINE OUTLIN TO OUTPUT RESULTS 

200 CALL OUTLIN 

210 END 
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SUBROUTINE SYM2A 


* THIS SUBROUTINE COMPUTES THE A MATRIX AND F VECTOR FOR MAIN 

* PROGRAM LIN2 AND SOLVES THE LINEAR SET OF EQUATIONS AU = F. 


SUBROUTINE SYM2A 

COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),COORD(100), 
ELEN.ICORRC100,2),NEL,NSNP 

DIMENSION AE(2,2), FS1E(2), FS2E(2), WKAREA(40600) 

ZERO OUT A MATRIX AND F VECTOR 

DO 210 IZ = 1,NSNP 
FT(IZ) = 0. 

DO 200 JZ = 1,NSNP 
A(IZ,JZ) = 0. 

CONTINUE 

CONTINUE 

ALPHA=0. 

PRINT*,'WHAT TYPE OF FORCING FUNCTION IS TO BE USED?’ 

PRINT*,’MIDPOINT = 1; 1/4 -3/4 APPROX = 2; CONSISTENT = 3' 

READ(6,*) NFF 

ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 
DO 375 IEL=1,NEL 


DETERMINE ELEMENTAL A MATRIX AND ELEMENTAL F VECTOR 


220 AE(1,1)=1./ELEN 

230 AE(1,2)=(-1. /ELEN) 

240 AE(2,1)=AE(1,2) 

250 AE(2,2)=AE(1,1) 

263 IF (NFF. EQ. 1) THEN 

264 FS1E(1)=3. *ELEN*(ALPHA+ELEN/2. ) 

265 FS1E(2)=FS1E(1) 

266 ELSEIF (NFF.EQ. 2) THEN 

267 FSlE(l)=3.*ELEN*(ALPHA+ELEN/4. ) 

268 FSlE(2)=3.*ELEN*(ALPHA+3.*ELEN/4. ) 

269 ELSE 

270 FS1E(l)-3. *ALPHA*ELEN+ELEN**2 

271 FS1E(2)=3. *ALPHA*ELEN+2. *(ELEN**2) 

272 ENDIF 


C DISTRIBUTE AE MATRICES AND FE VECTORS INTO SYSTEM A MATRIX AND F 
C VECTOR 


300 DO 370 11=1,2 

310 DO 350 JJ=1,2 

320 IN=ICORR(IEL,II) 

330 JN=IC0RR(IEL,JJ) 

340 A(IN,JN)=A(IN,JN) - AE(II,JJ) 

350 CONTINUE 

360 FT(IN)=FS1E(II) + FT(IN) 
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370 CONTINUE 

372 ALPHA=ALPHA + ELEN 

375 CONTINUE 

C IMPOSE KINEMATIC AND NATURAL BOUNDARY CONDITIONS 

376 A(l,l)=l. 

377 A(l,2)=0. 

378 FT(1)=0. 

379 FT(NSNP)=FT(NSNP)-12. 

380 M=1 

390 IDGT=3 
400 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE AU = F 

410 CALL LEQT2F(A,M,NSNP,IQ,FT,IDGTjWKAREA,IER) 

420 DO 440 NEW=1,NSNP 
430 U(NEVf)=FT(NEW) 

440 CONTINUE 
450 RETURN 
460 END 
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SUBROUTINE UX3EXT 


* 

ir 


c 

c 

c 

c 

c 


* 

★ 




* THIS SUBROUTINE COMPUTES THE VALUE OF U=X**3 AT THE SPECIFIED 

* NODAL POINTS FOR MAIN PROGRAM LIN2. 


* 

* 


100 SUBROUTINE UX3EXT 

110 COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),COORD(100), 
: ELEN,ICORR(100,2),NEL,NSNP 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) = COORD(NN)**3 

150 CONTINUE 
160 RETURN 
170 END 
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SOLUTION OF U" = 6X USING LUMPED MIDPOINT FORCING FUNCTION 


X-COORD 

U EXACT 

U FEM 

% DIFF 

0. 000 

0. 0000 

0.0000 

0. 0 

0. 200 

0. 0080 

0. 0120 

50. 0 

0.400 

0.0640 

0. 0720 

12. 5 

0.600 

0.2160 

0. 2280 

5. 6 

0. 800 

0. 5120 

0.5280 

3. 1 

1. 000 

1. 0000 

1. 0200 

2. 0 

1. 200 

1. 7280 

1. 7520 

1. 4 

1.400 

2. 7440 

2.7720 

1. 0 

1. 600 

4.0960 

4.1280 

0. 8 

1. 800 

5.8320 

5. 8680 

0. 6 

2. 000 

8.0000 

8.0400 

0.5 


* 
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SOLUTION OF U" = 6X USING 1/4 


X-COORD 

U EXACT 

U FEM 

0. 000 

0.0000 

0.0000 

0. 200 

0. 0080 

0. 0060 

0. 400 

0. 0640 

0. 0600 

0. 600 

0.2160 

0. 2100 

0. 800 

0.5120 

0. 5040 

1.000 

1. 0000 

0.9900 

1. 200 

1. 7280 

1. 7160 

1.400 

2.7440 

2.7300 

1.600 

4.0960 

4. 0800 

1. 800 

5.8320 

5.8140 

2.000 

8.0000 

7.9800 


3/4 LUMPED FORCING FUNCTION 
% DIFF 
0 . 0 
-25. 0 
-6. 3 
-2. 8 
- 1 . 6 
- 1 . 0 
-0. 7 
- 0.5 
-0.4 
-0. 3 
- 0 . 2 
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SOLUTION OF U" = 6X USING CONSISTENT FORCING FUNCTION 


X-COORD 

U EXACT 

U FEM 

% DIFF 

0. 000 

0. 0000 

0.0000 

0. 0 

0. 200 

0.0080 

0.0080 

0. 0 

0.400 

0.0640 

0. 0640 

0. 0 

0. 600 

0.2160 

0. 2160 

0. 0 

0. 800 

0.5120 

0. 5120 

0. 0 

1. 000 

1.0000 

1. 0000 

0. 0 

1. 200 

1. 7280 

1. 7280 

0. 0 

1.400 

2.7440 

2.7440 

0. 0 

1. 600 

4.0960 

4.0960 

0. 0 

1. 800 

5.8320 

5.8320 

0.0 

2.000 

8.0000 

8. 0000 

0. 0 
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200 CALL OUTLIN 
210 END 
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SUBROUTINE SYM3A 


C 

c 

c 

c 

c 

c 


* 

* 


* 

* 


* THIS SUBROUTINE COMPUTES THE A MATRIX AND F VECTOR FOR MAIN 

* PROGRAM LIN3 AND SOLVES THE LINEAR SET OF EQUATIONS AU = F. 

1 *--♦ —«_-» t ».->—i—» « _t> i i j l , i . x . . y - . ViV-A--A— t ,— iV - i t < - i V- i V -i*i—A--A—A**A—iV-A-iV-i*!—.ViV iVi l r-i 1 fA*‘i V , *A** l i*‘i l *‘ 


100 SUBROUTINE SYM3A 

110 COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),COORD(100), 
: ELEN,ICORR(100,2),NEL,NSNP 

120 DIMENSION AE(2,2), FS1E(2), FS2E(2), WKAREA(40600) 

C ZERO OUT A MATRIX AND F VECTOR 


140 DO 210 IZ = 1,NSNP 

150 FT(IZ) = 0. 

160 DO 200 JZ = 1,NSNP 

170 A(IZ,JZ) = 0. 

200 CONTINUE 

210 CONTINUE 

213 ALPHA=0. 

214 PRINT*,’WHAT TYPE OF FORCING FUNCTION IS TO BE USED?’ 

215 PRINT*,'MIDPOINT = 1; 1/4 -3/4 APPROX = 2; CONSISTENT = 3' 

216 READ(6,*) NFF 

C ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 


217 DO 375 IEL=1,NEL 

C DETERMINE ELEMENTAL A MATRIX AND ELEMENTAL F VECTOR 

220 AE(1,1)=1./ELEN 

230 AE(1,2)=(-1. /ELEN) 

240 AE(2,1)=AE(1,2) 

250 AE(2,2)=AE(1,1) 

263 IF (NFF. EQ. 1) THEN 

264 FS1E(1)=6. *ELEN*(ALPHA+ELEN/2. )**2 

265 FS1E(2)=FS1E(1) 

266 ELSEIF (NFF.EQ. 2) THEN 

267 FS1E(1)=6. *ELEN*(ALPHA+ELEN/4. )**2 

268 FS1E(2)=6. *ELEN*(ALPHA+3. *ELEN/4. )**2 

269 ELSE 

270 FS1E(1)=6. *(ALPHA**2)*ELEN+4.*ALPHA*(ELEN**2)+ELEN**3 

271 FS1E(2)=6. *(ALPHA**2)*ELEN+8.*ALPHA*(ELEN**2)+3. *ELEN**3 

272 ENDIF 

C DISTRIBUTE AE MATRICES AND FE VECTORS INTO SYSTEM A MATRIX AND F 
C VECTOR 


300 DO 370 11=1,2 

310 DO 350 JJ=1,2 

320 IN=ICORR(IEL.II) 

330 JN=ICORR(IEL,JJ) 

340 A(IN,JN)=A(IN,JN) - AE(II,JJ) 

350 CONTINUE 

360 FT(IN)=FS1E(II) + FT(IN) 
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370 CONTINUE 

372 ALPHA=ALPHA + ELEN 

375 CONTINUE 

C IMPOSE KINEMATIC AND NATURAL BOUNDARY CONDITIONS 

376 A(1,1)=1. 

377 A(1,2)=0. 

378 FT(1)=0. 

379 FT(NSNP)=FT(NSNP)-32. 

380 M=1 

390 IDGT=3 
400 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE AU = F 

410 CALL LEQT2F(A,M,NSNP,IQ,FT,IDGT,WKAREA,IER) 

420 DO 440 NEW=1,NSNP 
430 U(NEW)=FT(NEW) 

440 CONTINUE 
450 RETURN 
460 END 
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noonno 


* 

* 


SUBROUTINE UX4EXT »,* * 

’ * 


* THIS SUBROUTINE COMPUTES THE VALUE OF U=X**4 AT THE SPECIFIED * 

* NODAL POINTS FOR MAIN PROGRAM LIN3. * 

****************************************************************** 


100 SUBROUTINE UX4EXT 

110 COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),COORD(100), 
: ELEN,IC0RR(100,2),NEL,NSNP 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) = C00RD(NN)**4 

150 CONTINUE 
160 RETURN 
170 END 




SOLUTION OF U" = 12X**2 USING MIDPOINT LUMPED FORCING FUNCTION 


X-COORD 

U EXACT 

U FEM 

% DIFF 

0. 000 

0.0000 

0.0000 

0. 0 

0. 200 

0. 0016 

0. 0184 

1050. 3 

0.400 

0.0256 

0. 0608 

137. 5 

0. 600 

0.1296 

0.1848 

42. 6 

0. 800 

0.4096 

0.4864 

18. 8 

1. 000 

1. 0000 

1.1000 

10. 0 

1.200 

2.0736 

2.1984 

6. 0 

1.400 

3. 8416 

3. 9928 

3. 9 

1. 600 

6.5536 

6.7328 

2.7 

1.800 

10.4976 

10.7064 

2. 0 

2. 000 

16. 0000 

16.2400 

1.5 
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SOLUTION OF U" = 12X**2 USING 1/4 - 3/4 LUMPED FORCING FUNCTION 


X-COORD 

U EXACT 

U FEM 

% DIFF 

0. 000 

0.0000 

0.0000 

0. 0 

0. 200 

0.0016 

0.0046 

187. 8 

0.400 

0.0256 

0. 0296 

15. 7 

0. 600 

0.1296 

0.1326 

2. 3 

0. 800 

0.4096 

0.4096 

0. 0 

1.000 

1. 0000 

0.9950 

-0. 5 

1. 200 

2.0736 

2.0616 

-0. 6 

1.400 

3. 8416 

3.8206 

-0. 5 

1.600 

6.5536 

6.5216 

-0. 5 

1. 800 

10.4976 

10.4526 

-0.4 

2. 000 

16.0000 

15.9400 

-0.4 
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SOLUTION OF U" = 12X**2 USING CONSISTENT FORCING FUNCTION 


X-COORD 

U EXACT 

U FEM 

% DIFF 

0. 000 

0. 0000 

0. 0000 

0. 0 

0. 200 

0. 0016 

0.0016 

0. 1 

0. 400 

0. 0256 

0.0256 

0. 0 

0. 600 

0.1296 

0. 1296 

0. 0 

0. 800 

0.4096 

0.4096 

0. 0 

1. 000 

1. 0000 

1.0000 

0. 0 

1. 200 

2.0736 

2.0736 

0. 0 

1.400 

3.8416 

3.8416 

0. 0 

1. 600 

6.5536 

6.5536 

0. 0 

1. 800 

10.4976 

10.4976 

0.0 

2. 000 

16. 0000 

16.0000 

0 . 0 
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SUBROUTINE OUTLIN 


C 

c 

c 

c 

c 

c 

c 


★ 

* 


* THIS SUBROUTINE COMPUTES THE PERCENT ERROR BETWEEN THE EXACT 

* AND FEM VALUES OF U AND OUTPUTS THE APPROPRIATE INFORMATION 

* FOR MAIN PROGRAMS LIN|, LIN2, AND LIN3. 


* 

* 

* 

* 

* 


100 SUBROUTINE OUTLIN 

110 COMMON A(100,100),FT(100),U(100),UEXT(100),UDIF(100),COORD(100), 
:ELEN,ICORR(100,2),NEL,NSNP 

C COMPUTE PERCENT ERROR AT EACH NODAL POINT 


120 DO 150 IK=2,NSNP 

130 UDIF(IK)=100. *(U(IK)-UEXT(IK))/UEXT(IK) 

150 CONTINUE 

160 UDIF(1)=U(1)-UEXT(1) 

C OUTPUT RESULTS TO THE SCREEN AND TO A FILE 


170 WRITE(6,180) 

175 WRITE(30,180) 

180 FORMAT(/,IX,'X-COORD ' ,4X,’U EXACT’,4X,’U FEM',6X,’% DIFF’) 
190 WRITE(6,200) (COORD(NP).UEXT(NP),U(NP).UDIF(NP), NP=1,NSNP) 

195 WRITE(30,200) (COORD(NP),UEXT(NP),U(NP),UDIF(NP), NP=1,NSNP) 

200 FORMAT(/,2X,F5. 3,5X,F7. 4,3X,F7. 4,4X,F6. 1) 

230 RETURN 
240 END 
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APPENDIX C. LINEARIZATION VECTORS FOR CONSTANT 
LINEARIZATION TECHNIQUE 

For this analysis, h(u') in equation (3.6) is replaced by (u’f and the element 
linearization vector becomes 



(C.l) 


The three approximations of this integral outlined in III.B.l are now determined. 



Replacing h(u) in equation (3.7) with (u*) J yields 



(C.2) 



(C.3) 
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Linear Approximation 

Replacing h(u) with (u') J in equation (3.9.b) gives 
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APPENDIX D. LINEARIZATION MATRICES FOR THE CLASSICAL 

LINEARIZATION TECHNIQUE 

For this analysis, h{u') in equation (3.10) is replaced by u and the element 
linearization matrix integral becomes 


V- T » 
gg “ (■ 

J 0 


The two approximations of this integral outlined in III.B.2 are now determined. 


Midpoint Lumped Approximation 


Replacing h{u) in equation (3.11) with u' yields 


( u 0i+(Vt)i 


Linear Approximation 

Replacing h{u) with u in equation (3.12) gives 


J_ ± 

3 6 

i. j_ 

6 3 


(D. 2) 


f / * ( ^^ ) a*\/* \ f 

"Jo i|i_ij \ M \ /e J + (Uj+,) \4 


(D3) 


4 3 ( M j )i + ( w j+ l)i ( M j )i ( u j + l)i 

12 W)i + W+l)i ( w ‘). + 3 ( u j*+l)i 
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APPENDIX E. PROGRAM LISTINGS FOR CONSTANT 

LINEARIZATION 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

110 COMMON A(100,100),FS(100),FU(100),FT(100),U(100),U0LD(100), 

: UEXTC100),UDIF(100),C00RD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,IC0RR(100,2),NEL,NSNP,ITYPE 
115 C0NV=. 0001 

C READ IN PARAMETERS FROM DATA FILE 
130 READ(29,*) NEL,TLEN,COORD(1),ULBC,ITYPE,URBC 

C CALCULATE NUMBER OF NODAL POINTS 

135 NSNP=NEL+1 

C DETERMINE ELEMENT SIZE OF EQUAL LENGTHS 

137 ELEN=TLEN/FLOAT(NEL) 

C ESTABLISH LOCAL TO GLOBAL CORRESPONDENCE AND X COORDINATE OF 

C EACH NODE 

160 DO 169 IEIr=l,NEL 

162 ICORR(IEL,1)=IEL 

163 ICORR(IEL,2)=IEL+l 

165 COORD(IEL+l)=COORD(IEL) + ELEN 

169 CONTINUE 

C CALL SUBROUTINE NU2CAM TO CREATE A MATRIX AND F VECTOR 

170 CALL NU2CAM 

C CALL SUBROUTINE NU2CAI TO PERFORM SOLUTION ITERATION 

180 CALL NU2CAI(IET) 


PROGRAM NU2CAN * 

THIS PROGRAM SOLVES THE NONLINEAR SECOND ORDER DIFFERENTIAL * 
EQUATION: * 

U" - U**2 = 6 - 9X**4; UEXACT=3X**2 WITH VARIABLE DOMAIN * 
TREATING THE U**2 TERM AS AN EXCITATION AND TAKING IT TO THE * 
RIGHT SIDE OF THE EQUATION. * 

THE USER SELECTS: * 

1) NUMBER OF ELEMENTS * 

2) SIZE OF DOMAIN * 

3) X AND U(X) AT THE LEFT BOUNDARY * 

4) U(X) OR U'(X) AT THE RIGHT BOUNDARY * 

5) ITERATION STRATEGY FOR DETERMINING U* * 

6) APPROXIMATION TECHNIQUE FOR THE EXCITATION INTEGRAL * 
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C CALL SUBROUTINE U2EXTA TO COMPUTE EXACT SOLUTION U=3X**2 
190 CALL U2EXTA 

C CALL SUBROUTINE OUTPUT TO PRINT OUT DATA, COMPUTATIONAL EFFICIENCY 

200 CALL OUTPUT(CPUSTAR,IET) 

210 END 
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c 

c 

c 

c 

c 

c 



* THIS SUBROUTINE COMPUTES THE A MATRIX ANT F VECTOR FOR MAIN * 

* PROGRAM NU2CAN. * 


100 SUBROUTINE NU2CAM 

110 COMMON A(100,100),FS(100),FU(100),FT(100),U(100),UOLD(100), 

: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME.ULBC.URBC, 

: TLEN,ICORR(100,2),NEL,NSNP,ITYPE 
120 DIMENSION AE(2,2), FS1E(2), FS2E(2) 

122 IF (TLEN. LE. 1. 0) THEN 

123 PRINT*,'CHOOSE BOUNDARY FOR INITIAL GUESS:' 

124 PRINT*,'1 = LEFT ESSENTIAL BOUNDARY CONDITION' 

125 PRINT*,'2 = RIGHT ESSENTIAL BOUNDARY CONDITION' 

126 PRINT*,'3 = AVERAGE OF THE TWO ESSENTIAL BOUNDARY CONDITIONS’ 

127 READ(6,*) INITGS 

128 ELSE 

129 CONTINUE 

130 ENDIF 

140 DO 210 IZ = 1,NSNP 


C ZERO OUT STEADY FORCE VECTOR 


150 FS(IZ) = 0. 

C DETERMINE INITIAL VALUE OF USTAR TO BEGIN THE ITERATION PROCESS 


157 IF (INITGS. EQ.1) THEN 

158 U(IZ)=ULBC 

159 UOLD(IZ)=ULBC 

160 ELSEIF (INITGS. EQ. 2) THEN 

161 U(IZ)=URBC 

162 UOLD(IZ)=URBC 

163 ELSEIF (INITGS. EQ.3) THEN 

164 U(IZ)=(ULBC+URBC)/2. 

165 UOLD(IZ)=U(IZ) 

166 ELSE 

167 U(IZ)=SQRT(ABS(9.*C00RD(IZ)**4 - 6.)) 

168 UOLD(IZ)=U(IZ) 

169 ENDIF 


C ZERO OUT A MATRIX 

170 DO 200 JZ = 1,NSNP 

180 A(IZ,JZ) = 0. 

200 CONTINUE 

210 CONTINUE 

C ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 


213 ALPHA=0. 

215 DO 375 IEL=1,NEL 
220 AE(1,1)=1. /ELEN 
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230 

AE(1,2)=(-1. /ELEN) 


240 

AE(2,1)=AE(1,2) 


250 

AE(2,2)=AE(1,1) 


260 

FS1E(1)=3.*ELEN 


270 

FS1E(2)=FS1E(1) 


272 

Fl=(ALPHA**4)*ELEN/2. 


274 

F2=2. *(ALPHA**3)*(ELEN**2)/3. 


276 

F3=(ALPHA**2)*(ELEN**3)/2. 


278 

F4=ALPHA*(ELEN**4)/5. 


280 

F5=(ELEN**5)/30. 


287 

FS2E(l)=(-9.)*(F1 + F2 + F3 + F4 

+ F5) 

290 

FS2E(2)=(-9)*(F1 + 2.*F2 + 3. *F3 

+ 4. *F4 

300 

DO 370 11=1,2 


310 

DO 350 JJ=1,2 


320 

IN=ICORR(IEL,II) 


330 

JN=ICORR(IEL,JJ) 


340 

A(IN,JN)=A(IN,JN) - AE(II,JJ) 

350 

CONTINUE 


360 

FS(IN)=FS1E(II) + FS2ECII) + 

FS(IN) 

370 

CONTINUE 

372 

ALPHA=ALPHA + ELEN 


375 

CONTINUE 


420 

RETURN 


430 

END 



5.*F5) 
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SUBROUTINE NU2CAI 


C 

c 

c 

c 

c 

c 


★ 

* 


•kitickir/tirkir/titiririritirkirkirlrirkicirlrk'ic 

* 

* 


* THIS SUBROUTINE PERFORMS THE ITERATIVE SOLUTION PROCESS FOR 

* MAIN PROGRAM NU2CAN. 


100 SUBROUTINE NU2CAI(IET) 

102 COMMON A(100,100),FS(100),FU(100),FT(100),U(100),UOLD(100), 

: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORRv100,2),NEL,NSNP,ITYPE 

104 DIMENSION WKAREAC40600), DIF(IOO), FUE(2), USTAR(IOO) 

C SELECT METHOD OF DETERMINING USTAR 

105 PRINT*,'SELECT METHOD OF U* DETERMINATION. ' 

106 PRINT*,'1: U* = U* 

107 PRINT*,'2: U* = AVERAGE OF LAST TWO COMPUTED VALUES OF U 
109 READ(6,*) METHU 

C SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR 

116 PRINT*,'SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR.' 

117 PRINT*,'1: MIDPOINT APPROXIMATION’ 

118 PRINT*,'2: 1/4 - 3/4 APPROXIMATION* 

119 PRINT*,'3: LINEAR' 

120 READ(6,*) METHFU 

C CALL SUBROUTINE SETIME TO BEGIN TIMING ITERATION PROCESS 


121 CALL SETIME 

C BEGIN ITERATION PROCESS 

122 DO 450 ITER=1,200 

C RESET VALUE OF UNSTEADY F VECTOR TO ZERO 

123 DO 138 IU=1,NSNP 

124 FU(IU)=0. 

C DETERMINE VALUE OF U* AT EACH NODE 

125 IF (METHU. EQ. 1) THEN 

126 USTAR(IU)=U(IU) 

127 ELSEIF (METHU. EQ. 2) THEN 

128 USTARCIU)=(U(IU)+UOLD(IU))/2. 

132 ENDIF 

138 CONTINUE 

C DETERMINE UNSTEADY FORCE VECTOR 

140 DO 210 IEL=1,NEL 

145 IF (METHFU. EQ. 1) THEN 

146 FUE(1)=(ELEN/2.)*((USTAR(IEL)+USTAR(IEL+l))/2.)**2 

147 FUE(2)=FUE(1) 
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149 ELSEIF (METHFU. EQ.2) THEN 

151 FUE(l)=(ELEN/2. )*(3. *USTAR(IEL)/4. + USTAR(IEL+l)/4. )**2 

152 FUE(2)=(ELEN/2. )*(USTAR(IEL)/4. + 3. *USTAR(IEL+l)/4. )**2 

153 ELSE 

154 FUE(l)=ELEN*(USTAR(IEL)**2/4. +USTARCIEL)*USTAR(IEL+l)/6. 

: + USTAR(IEL+1)**2/12. ) 

155 FUE(2)=ELEN*(USTAR(IEL)**2/12. + USTAR(IEL)*USTAR(IEL+1) 

: /6. + USTAR(IEL+1)**2/4.) 

156 ENDIF 

170 DO 200 11=1,2 

180 IN=ICORR(IEL,II) 

190 FU(IN)=FUE(II) + FU(IN) 

200 CONTINUE 

210 CONTINUE 

C DETERMINE TOTAL FORCE VECTOR 

220 DO 240 NP=1,NSNP 

230 FT(NP)=FS(NP)+FU(NP) 

235 UOLD(NP)=U(NP) 

240 CONTINUE 

C IMPOSE BOUNDARY CONDITIONS 

241 A(l,l)-1. 

242 A(l,2)=0. 

243 FT(1)=ULBC 

244 IF (ITYPE.EQ. 1) THEN 

245 A(NSNP,NSNP-1)=0. 

246 A(NSNP,NSNP)=1. 

247 FT(NSNP)=URBC 

248 ELSE 

249 FT(NSNP)=FT(NSNP)-URBC+TLEN 

250 ENDIF 

255 M=1 

260 IDGT=3 

270 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE SET OF LINEAR ALGEBRAIC EQUATIONS 

280 CALL LEQT2F(A,M,NSNP,IQ,FT,IDGT,WKAREA,IER) 

290 DO 310 NEW=1,NSNP 

300 U(NEW)=FT(NEW) 

C WRITEC*,*) 'UNEV=',U(NEW) 

C TEST FOR CONVERGENCE 

306 DIF(NEW)=ABS(U(NEW)-UOLD(NEW)) 

310 CONTINUE 

320 DIFMAX=DIF(1) 

325 NMAX=1 

330 DO 390 1J=1,NEL 

340 IF (DIF(IJ+1).GE.DIF(IJ)) THEN 

350 DIFMAX=DIF(IJ+1) 

355 NMAX=IJ+1 

360 ELSE 
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370 

CONTINUE 


380 

ENDIF 


390 

CONTINUE 

405 

IF (ABS(DIFMAX/U(NMAX)). LT. CONV) THEN 


410 

GO TO 451 


420 

ELSE 


430 

CONTINUE 


440 

ENDIF 


450 

CONTINUE 


451 

CONTINUE 


C 

CALL SUBROUTINE GETIME TO OBTAIN CPU TIME FOR ITERATION PROCESS 

452 

CALL GETIME(IET) 


C 

OUTPUT HEADER INFORMATION 


462 

WRITEC6,464) 


463 

WRITEC 30,464) 

FORMATCIX,'EQUATION: U" - U**2 - 6 - 9X**4') 


464 


465 

IF (ITYPE.EQ. 1) THEN 


466 

WRITEC6,468) COORD(1),ULBC,COORD(NSNP),URBC 


467 

WRITE(30,468) COORD(1),ULBC,COORD(NSNP),URBC 

FORMATCIX,'B.C.: U(',F2. 0,')=',F2. 0,'; U(',F3.0, , )= 


468 

',F4.0,/) 

469 

ELSE 


470 

WRITE(6,472) COORD(1),ULBC,COORD(NSNP),URBC 


471 

WRITEC 30,472) COORD(1),ULBC,COORD(NSNP),URBC 

FORMATCIX,'B.C.: UC',F2. 0,')=’,F2. 0,'; DU/DXC',F3. 0, 

’)=’,F4.0,/) 

472 

473 

ENDIF 


475 

IF CMETHU.EQ. 1) THEN 


476 

WRITEC6.478) 


477 

WRITEC30.478) 

FORMATCIX,'ITERATION METHOD: U*=U',/) 


478 


479 

ELSE 


480 

WRITEC6.482) 


481 

WRITEC30,482) 

FORMATCIX,'ITERATION METHOD: U*=CU+U0LD)/2’,/) 


482 


483 

ENDIF 


488 

IF CMETHFU.EQ. 1) THEN 


489 

WRITEC6.491) 


490 

WRITEC30,491) 

FORMATCIX,’METHOD OF EXCITATION INTEGRAL EVALUATION: 


491 

MIDPOINT',/) 

492 

ELSEIF CMETHFU.EQ. 2) THEN 


493 

WRITEC6,495) 


494 

WRITEC30,495) 

FORMATCIX,’METHOD OF EXCITATION INTEGRAL EVALUATION: 


495 

1/4-3/4',/) 

496 

ELSE 


497 

WRITEC6.499) 


498 

WRITEC30,499) 


499 

FORMATCIX,'METHOD OF EXCITATION INTEGRAL EVALUATION: 

LINEAR',/) 

500 

ENDIF 


502 

IF CITER.GE. 200) THEN 


503 

WRITEC6.505) 


504 

WRITEC30,505) 


505 

FORMATCIX, CONVERGENCE NOT OBTAINED AFTER 200 ITERATIONS. ') 

506 

ELSEIF CABSCUCNMAX)).GT. CIO.**20). OR.ABSCUCNSNP-1)). GT. 
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:(10.**20)) THEN 

507 WRITE(6,509) 

508 WRITE(30,509) 

509 FORMAT(IX,'SOLUTION PROCESS DIVERGES.') 

510 ELSE 

511 WRITE(6,520) ITER,NEL 

515 WRITE(30,520) ITER,NEL 

520 FORMAT(IX,'CONVERGENCE OBTAINED AFTER ',13,' ITERATIONS USING ', 

: 13,’ ELEMENTS. ',/) 

525 ENDIF 
530 RETURN 
540 END 
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c 

c 

c 

c 

c 

c 



100 SUBROUTINE U2EXTA 

110 COMMON A(100,100),FS(100),FU(100),FT(100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),NEL,NSNP,ITYPE 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) = 3. *C00RD(NN)**2 

150 CONTINUE 
160 RETURN 
170 END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


* 

* 

* 

★ 

* 

* 

* 

* 

* 

* 

* 

* 


PROGRAM NU2CBN * 

THIS PROGRAM SOLVES THE NONLINEAR SECOND ORDER DIFFERENTIAL * 

EQUATION: * 

U" + U**2 = 60X + 100X**6; UEXACT=10X**3 WITH VARIABLE * 
DOMAIN, TREATING THE U**2 TERM AS AN EXCITATION AND TAKING IT * 
TO THE RIGHT SIDE OF THE EQUATION. THE USER SELECTS: * 

1) NUMBER OF ELEMENTS * 

SIZE OF DOMAIN * 

X AND U(X) AT THE LEFT BOUNDARY * 

U(X) OR U'(X) AT THE RIGHT BOUNDARY * 

ITERATION STRATEGY FOR DETERMINING U* * 

APPROXIMATION TECHNIQUE FOR THE EXCITATION INTEGRAL * 


2 ) 

3) 

4) 

5) 

6 ) 


110 COMMON A(100,100),FS(100),FU(100),FT(100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),NEL,NSNP,ITYPE 

115 CONV=. 0001 


C READ IN PARAMETERS FROM DATA FILE 

130 READ(29,*) NEL,TLEN,COORD(1),ULBC,ITYPE,URBC 

C CALCULATE NUMBER OF NODAL POINTS 


135 NSNP=NEL+1 

C DETERMINE ELEMENT SIZE OF EQUAL LENGTHS 

137 ELEN=TLEN/FLOAT(NEL) 

C ESTABLISH LOCAL TO GLOBAL CORRESPONDENCE AND X COORDINATE OF 
C EACH NODE 


160 DO 169 IEL=1,NEL 

162 ICORR(IEL,1)=IEL 

163 IC0RR(IEL,2)=IEL+1 

164 C00RD(IEL+1)=C00RD(IEL)+ELEN 

169 CONTINUE 

C CALL SUBROUTINE NU2CBM TO CREATE A MATRIX AND F VECTOR 

170 CALL NU2CBM 

C CALL SUBROUTINE NU2CBI TO PERFORM SOLUTION ITERATION 

180 CALL NU2CBI(IET) 

C CALL SUBROUTINE U2EXTB TO COMPUTE EXACT SOLUTION U=10X**3 

190 CALL U2EXTB 

C CALL SUBROUTINE OUTPUT TO PRINT OUT DATA, COMPUTATIONAL EFFICIENCY 
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200 CALL OUTPUT(CPUSTAR,IET) 
210 END 





* 


100 SUBROUTINE NU2CBM 

110 COMMON A(100,100),FS(100),FU(100),FT(100),U(100),UOLD(100), 

: UEXTC100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),NEL.NSNP,ITYPE 
120 DIMENSION AE(2,2), FS1E(2), FS2E(2) 

122 IF (TLEN. LT. 1.0) THEN 

123 PRINT*,'CHOOSE BOUNDARY FOR INITIAL GUESS:' 

124 PRINT*,'1 = LEFT ESSENTIAL BOUNDARY CONDITION' 

125 PRINT*,'2 = RIGHT ESSENTIAL BOUNDARY CONDITION' 

126 PRINT*,'3 = AVERAGE OF THE TWO ESSENTIAL BOUNDARY CONDITIONS' 

127 READ(6,*) INITGS 

128 ELSE 

129 CONTINUE 

130 ENDIF 

140 DO 210 IZ = l.NSNP 
C ZERO OUT STEADY FORCE VECTOR 

150 FS(IZ) = 0. 

C DETERMINE INITIAL VALUE OF USTAR TO BEGIN THE ITERATION PROCESS 

157 IF (INITGS. EQ. 1) THEN 

158 U(IZ)=ULBC 

159 UOLD(IZ)=ULBC 

160 ELSEIF (INITGS. EQ. 2) THEN 

161 U(IZ)=URBC 

162 UOLD(IZ)=URBC 

163 ELSEIF (INITGS. EQ. 3) THEN 

164 U(IZ)=(ULBC+URBC)/2. 

165 UOLD(IZ)=U(IZ) 

166 ELSE 

167 U(IZ)=SQRT(60.*C00RD(IZ) + 100.*COORD(IZ)**6) 

168 UOLD(IZ)=U(IZ) 

169 ENDIF 

C ZERO OUT A MATRIX 

170 DO 200 JZ = l.NSNP 

175 A(IZ,JZ) = 0. 

200 CONTINUE 

210 CONTINUE 

C ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 


SUBROUTINE NU2CBM 


* THIS SUBROUTINE COMPUTES THE A MATRIX AND F VECTOR FOR MAIN 

* PROGRAM NU2CBN. 


213 ALPHA=0. 

215 DO 375 IEL=1,NEL 
220 AE(1,1)=1./ELEN 




230 AE(1,2)=(-1. /ELEN) 

240 AE(2,1)=AE(1,2) 

250 AE(2,2)=AE(1,1) 

260 FS1E(1)=30.*ALPHA*ELEN + 10.*ELEN**2 

270 FS1E(2)=30.*ALPHA*ELEN + 20.*ELEN**2 

272 F1=50.*(ALPHA**6)*ELEN 

274 F2=100.*( ALPHA**5 )*(ELEN**2) 

276 F3=125.*(ALPHA**4)*(ELEN**3) 

278 F4=100.*(ALPHA**3)*(ELEN**4) 

280 F5=50.*(ALPHA**2)*(ELEN**5) 

281 F6=100.*ALPHA*(ELEN**6)/7. 

282 F7=25.*(ELEN**7)/16. 

287 FS2E(1)=F1 + F2 + F3 + F4 + F5 + F6 + F7 

290 FS2E(2)=F1 + 2.*F2 + 3.*F3 + 4. *F4 + 5.*F5 + 6.*F6 + 7*F7 

300 DO 370 11=1,2 

310 DO 350 JJ=1,2 

320 IN=ICORR(IEL,II) 

330 JN=ICORR(IEL,JJ) 

340 A(IN,JN)=A(IN,JN) - AE(II.JJ) 

350 CONTINUE 

360 FS(IN)=FS1E(II) + FS2E(II) + FS(IN) 

370 CONTINUE 

372 ALPHA=ALPHA + ELEN 

375 CONTINUE 

420 RETURN 

430 END 
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SUBROUTINE NU2CBI 


C 

c 

c 

c 

c 

c 


* 

* 


'•t 


* 

* 


* THIS SUBROUTINE PERFORMS THE ITERATIVE SOLUTION PROCESS FOR 

* MAIN PROGRAM NU2CBN. 


100 SUBROUTINE NU2CBI(IET) 

102 COMMON A(100,100),FS(100),FU(100),FT(100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),NEL,NSNP,ITYPE 
104 DIMENSION WKAREA(40600), DIF(IOO), FUE(2), USTAR(IOO) 


C SELECT METHOD OF DETERMINING USTAR 


105 PRINT*,'SELECT METHOD OF U* DETERMINATION. ' 

106 PRINT*,'1: U* = U' 

107 PRINT*,'2: U* = AVERAGE OF LAST TWO COMPUTED VALUES OF U' 

109 READ(6,*) METHU 


C SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR 


116 

117 

118 

119 

120 


PRINT*,'SELECT METHOD OF DETERMINING UNSTEADY FORCE 
PRINT*,'1: MIDPOINT APPROXIMATION' 

PRINT*,'2: 1/4 - 3/4 APPROXIMATION’ 

PRINT*,'3: CONSISTENT' 

READ(6,*) METHFU 


VECTOR. ' 


C CALL SUBROUTINE SETIME TO BEGIN TIMING ITERATION PROCESS 


121 CALL SETIME 


C BEGIN ITERATION PROCESS 


122 DO 450 ITER=1,200 

C RESET VALUE OF UNSTEADY F VECTOR TO ZERO 


123 DO 138 IU=1,NSNP 

124 FU(IU)=0. 

C DETERMINE VALUE OF U* AT EACH NODE 

125 IF (METHU. EQ. 1) THEN 

126 USTARCIU)=U(IU) 

127 ELSE 

128 USTARCIU)=(U(IU)+UOLD(IU))/2. 

132 ENDIF 

138 CONTINUE 

C DETERMINE UNSTEADY FORCE VECTOR 

140 DO 210 IEL=1,NEL 

145 IF (METHFU. EQ. 1) THEN 

146 FUE(1)=(ELEN/2. )*((USTARCIEL)+USTAR(IEL+1))/2. )**2 

147 FUE(2)=FUE(1) 
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149 ELSEIF (METHFU.EQ. 2) THEN 

151 FUE(1)=(ELEN/2. )*(3. *USTAR(IEL)/4. + USTAR(IEL+l)/4. )**2 

152 FUE(2)=(ELEN/2. )*(USTAR(IEL)/4. + 3. *USTAR(IEL+l)/4. )**2 

153 ELSE 

154 FUE(l)=ELEN*(USTAR(IEL)**2/4. +USTAR(IEL)*USTAR(IEL+l)/6. 

: + USTAR(IEL+1)**2/12. ) 

155 FUE(2)=ELEN*(USTAR(IEL)**2/12. + USTAR(IEL)*USTAR(IEL+l)/6. 

: + USTAR(IEL+l)**2/4. ) 

156 ENDIF 

170 DO 200 11=1,2 

180 IN=ICORR(IEL,II) 

190 FU( IN)=FUE(II) + FU(IN) 

200 CONTINUE 

210 CONTINUE 

C DETERMINE TOTAL FORCE VECTOR 

220 DO 240 NP=1,NSNP 

230 FT(NP)=FS(NP)-FU(NP) 

235 UOLD(NP)=U(NP) 

240 CONTINUE 

C IMPOSE BOUNDARY CONDITIONS 

241 A(1,1)=1. 

242 A(l,2)=0. 

243 FT(1)=ULBC 

244 IF (ITYPE.EQ. 1) THEN 

245 A(NSNP,NSNP-1)=0. 

246 A(NSNP,NSNP)=1. 

247 FT(NSNP)=URBC 

248 ELSE 

249 FT(NSNP)=FT(NSNP)-URBC 

250 ENDIF 

255 M=1 

260 IDGT=3 

270 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE SET OF LINEAR ALGEBRAIC EQUATIONS 

280 CALL LEQT2F(A,M,NSNP,IQ,FT,IDGT,WKAREA,IER) 

290 DO 310 NEW=1,NSNP 

300 U(NEW)=FT(NEW) 

C WRITE(*,*) 1 UNEW=',U(NEV) 

C TEST FOR CONVERGENCE 

306 DIF(NEW)=ABS(U(NEW)-UOLD(NEW)) 

310 CONTINUE 

320 DIFMAX=DIF(1) 

325 NMAX=1 

330 DO 390 IJ=1,NEL 

340 IF (DIF(IJ+1). GE. DIF(IJ)) THEN 

350 DIFMAX=DIF(IJ+1) 

355 NMAX=IJ+1 

360 ELSE 
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370 CONTINUE 

380 ENDIF 

390 CONTINUE 

405 IF (ABS(DIFMAX/U(NMAX)). LT. CONV) THEN 

410 GO TO 460 

420 ELSE 

430 CONTINUE 

440 ENDIF 

450 CONTINUE 

460 CONTINUE 

C CALL SUBROUTINE GETIME TO OBTAIN CPU TIME FOR ITERATION PROCESS 

461 CALL GETIME(IET) 

C OUTPUT HEADER INFORMATION 

462 WRITE(6,464) 

463 WRITE(30,464) 

464 FORMATCIX,'EQUATION: U" + U**2 *= 60X + 100X**6’) 

465 IF (ITYPE.EQ. 1) THEN 

466 WRITE(6,468) COORD(1),ULBC.COORD(NSNP),URBC 

467 WRITE(30,468) COORD(1),ULBC,COORD(NSNP),URBC 

468 FORMAT(IX,’B. C.: U(',F2.0,')=',F3.0,'; U(',F2. 0,’)=’,F4.0,/) 

469 ELSE 

470 WRITE(6,472) COORD(1),ULBC,COORD(NSNP),URBC 

471 WRITE(30,472) COORD(1),ULBC,COORD(NSNP),URBC 

472 FORMATCIX,'B. C.: UC',F2. 0,')=',F3. 0,’; DU/DXC’,F2. 0,')=’,F4. 0,/) 

473 ENDIF 

475 IF CMETHU.EQ. 1) THEN 

476 WRITEC6.478) 

477 WRITEC 30,478) 

478 FORMATCIX,'ITERATION METHOD: U*=U',/) 

479 ELSE 

480 WRITEC6.482) 

481 WRITEC30.482) 

482 FORMATC IX,'ITERATION METHOD: U*KU+U0LD)/2' ,/) 

487 ENDIF 

488 IF CMETHFU.EQ. 1) THEN 

489 WRITEC6,491) 

490 WRITEC 30,491) 

491 FORMATCIX,'METHOD OF EXCITATION INTEGRAL EVALUATION: MIDPOINT’,/) 

492 ELSEIF CMETHFU. EQ. 2) THEN 

493 WRITEC6.495) 

494 WRITEC30.495) 

495 FORMATCIX,'METHOD OF EXCITATION INTEGRAL EVALUATION: 1/4-3/4',/) 

496 ELSE 

497 WRITEC6.499) 

498 WRITEC30,499) 

499 FORMATCIX,'METHOD OF EXCITATION INTEGRAL EVALUATION: LINEAR',/) 

500 ENDIF 

502 IF (ITER. GE.200) THEN 

503 WRITEC6.505) 

504 WRITEC30,505) 

505 FORMATCIX,'CONVERGENCE NOT OBTAINED AFTER 200 ITERATIONS.’) 

506 ELSEIF CABSCUCNMAX)).GT.C10. **20).OR.ABSCUCNSNP-1)). GT. 
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:(10.**20)) THEN 

507 WRITE(6,509) 

508 WRITE(30,509) 

509 FORMATCIX,'SOLUTION PROCESS DIVERGES.') 

510 ELSE 

511 WRITE(6,520) ITER.NEL 

515 WRITE(30,520) ITER.NEL 

520 FORMATCIX,'CONVERGENCE OBTAINED AFTER ',13,' ITERATIONS USING ', 

:13,' ELEMENTS. ',/) 

525 ENDIF 
530 RETURN 
540 END 
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SUBROUTINE U2EXTB 


C 

c 

c 

c 

c 

c 


* 

* 


* 

* 


* THIS SUBROUTINE COMPUTES THE EXACT SOLUTION, U=10X**3, 

* FOR MAIN PROGRAM NU2CBN AT THE SPECIFIED NODAL POINTS. 


100 SUBROUTINE U2EXTB 

110 COMMON A(100,100),FS(100),FU(100),FT(100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),NEL,NSNP,ITYPE 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) * 10. *C00RD(NN)**3 

150 CONTINUE 
160 RETURN 
170 END 
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C * SUBROUTINE OUTPUT 

C * THIS SUBROUTINE COMPUTES THE PER CENT ERROR BETWEEN THE EXACi 
C * AND FEM SOLUTIONS, CPU* FOR THE ITERATION PROCESS, AND PRINTS 

C * OUT ALL DATA IN TABULAR FORM FOR PROGRAMS NU2CAN AND NU2CBN 

100 SUBROUTINE OUTPUT(CPUSTAR.IET) 

110 COMMON A(100,100),FS(100),FU(100),FT(100),U(500),UOLD(100), 

:UEXT(100),UDIF(100),COORDC100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),NEL,NSNP,ITYPE 
115 SUMDIF=0. 

C CALCULATE PER CENT ERROR AT EACH NODE AND SUM THE ABSOLUTE VALUE 
C OF ALL THE ERRORS 

120 DO 150 IK=2,NSNP 

130 UDIF(IK)=100. *(U(IK)-UEXT(IK))/UEXT(IK) 

140 SUMDIF=SUMDIF + ABS(UDIF(IK)) 

150 CONTINUE 

160 UDIF(1)=U(1)-UEXT(1) 

C COMPUTE THE ELAPSED TIME OF THE ITERATION PROCESS 

164 ELTIME=IET*. 000026 

165 WRITE(6,169) ELTIME 

166 WRITEC 30,169)ELTIME 

169 FORMAT(IX,'ELAPSED TIME FOR THE ITERATION PROCESS IS ’,F9.4, 

: ' SECONDS. ') 

C OUTPUT DATA IN TABULAR FORMAT 

170 WRITE(6,180) 

175 WRITE(30,180) 

180 FORMAT(/,IX,'X-COORD',3X,'U EXACT’,3X,’U FEM',4X,'% DIFF’) 

190 WRITE(6,200) (COORD(NP).UEXT(NP),U(NP),UDIF(NP), NP=1,NSNP) 

195 WRITE(30,200) (COORD(NP),UEXT(NP),U(NP).UDIF(NP), NP=1,NSNP) 

200 FORMAT(/,2X,F5. 3,5X,F7. 4,3X,F7. 4,4X,F5.1) 

C CALCULATE CPU* FOR THE ITERATION PROCESS 

205 CPUSTAR=ELTIME*SUMDIF/NSNP 
210 WRITE(6,220) CPUSTAR 

215 WRITEC30,220) CPUSTAR 

220 F0RMAT(/,1X, CPU* FOR THE ITERATION PROCESS IS ’,F9.4,' SECONDS.’) 
230 RETURN 

240 END 
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APPENDIX F. PROGRAM LISTINGS FOR CLASSICAL 

LINEARIZATION 

c *1 

c * 

c * 

c * 

c * 

c * 

c * 

c * 

c * 

c * 

c * 

c * 

c * 

c *• 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),U0LD(100), 
: UEXT(100),UDIF(100),C00RD(100),ELEN,C0NV,ELTIME,ULBC,URBC, 

: TLEN,IC0RR(100,2),ITYPE,NEL,NSNP 
115 C0NV=. 0001 

C READ IN PARAMETERS FROM DATA FILE 
130 READ(29,*) NEL,TLEN,COORD(1),ULBC,ITYPE,URBC 

C CALCULATE NUMBER OF NODAL POINTS 
135 NSNP=NEL+1 

C DETERMINE ELEMENT SIZE OF EQUAL LENGTHS 

137 ELEN=TLEN/FLOAT(NEL) 

C ESTABLISH LOCAL TO GLOBAL CORRESPONDENCE AND X COORDINATE OF 
C EACH NODE 

160 DO 169 IEL=1,NEL 

162 ICORR(IEL,1)=IEL 

163 ICORR(IEL,2)=IEL+l 

164 COORD(IEL+l)=COORD(IEL)+ELEN 

169 CONTINUE 

C CALL SUBROUTINE NU2CAM TO CREATE A MATRIX AND F VECTOR 

170 CALL NU2KAM 

C CALL SUBROUTINE NU2CAI TO PERFORM SOLUTION ITERATION 

180 CALL NU2KAI(IET) 


PROGRAM NU2KA * 

THIS PROGRAM SOLVES THE NONLINEAR SECOND ORDER DIFFERENTIAL * 
EQUATION: * 

U" - U**2 = 6 - 9X**4; UEXACT=3X**2 WITH VARIABLE DOMAIN * 

BY LINEARIZING THE U**2 TERM AS USTAR*U AND KEEPING IT ON THE * 

LEFT SIDE OF THE EQUATION. THE USER SELECTS: * 

1) NUMBER OF ELEMENTS * 

2) SIZE OF DOMAIN * 

3) X AND U(X) AT THE LEFT BOUNDARY * 

4) U(X) OR U'(X) AT THE RIGHT BOUNDARY * 

5) ITERATION STRATEGY FOR DETERMINING U* * 

6) APPROXIMATION TECHNIQUE FOR THE EXCITATION INTEGRAL * 
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C CALL SUBROUTINE U2EXTA TO COMPUTE EXACT SOLUTION U=3X**2 
190 CALL CLEXTA 

C CALL SUBROUTINE OUTPUT TO PRINT OUT DATA, COMPUTATIONAL EFFICIENCY 

200 CALL CLOTPT(CPUSTAR,IET) 

210 END 



c 

c 

c 

c 

c 

c 


* SUBROUTINE NU2KAM * 

* * 

* THIS SUBROUTINE COMPUTES THE A MATRIX AND F VECTOR FOR MAIN * 

* PROGRAM NU2KA. * 


100 SUBROUTINE NU2KAM 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
:UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),ITYPE,NEL,NSNP 
120 DIMENSION AE(2,2), FS1E(2), FS2E(2) 

122 IF (TLEN. LE. 1.0) THEN 

123 PRINT*,'CHOOSE BOUNDARY FOR INITIAL GUESS:’ 

124 PRINT*,'1 = LEFT ESSENTIAL BOUNDARY CONDITION’ 

125 PRINT*,'2 = RIGHT ESSENTIAL BOUNDARY CONDITION' 

126 PRINT*,'3 = AVERAGE OF THE TWO ESSENTIAL BOUNDARY CONDITIONS' 

127 READ(6,*) INITGS 

128 ELSE 

129 CONTINUE 

130 ENDIF 

140 DO 210 IZ = 1,'NSNP 

C ZERO OUT STEADY FORCE VECTOR 
150 FS(IZ) = 0. 

C DETERMINE INITIAL VALUE OF USTAR TO BEGIN THE ITERATION PROCESS 

157 IF (INITGS. EQ. 1) THEN 

158 U(IZ)=ULBC 

159 UOLD(IZ)=U(IZ) 

160 ELSEIF (INITGS. EQ. 2) THEN 

161 U(IZ)=URBC 

162 UOLD(IZ)=U(IZ) 

163 ELSEIF (INITGS. EQ. 3) THEN 

164 U(IZ)=(ULBC+URBC)/2. 

165 UOLD(IZ)=U(IZ) 

166 ELSE 

167 U(IZ) = SQRT(ABS(9.*COORD(IZ)**4 - 6.)) 

168 UOLD(IZ) = U(IZ) 

169 ENDIF 

C ZERO OUT ALL MATRICES 

170 DO 200 JZ = 1,NSNP 

171 A(IZ,JZ) = 0. 

175 B(IZ,JZ) = 0. 

176 C(IZ,JZ) = 0. 

200 CONTINUE 

210 CONTINUE 

C ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 

213 ALPHA=0. 

215 DO 375 IEL=1,NEL 
220 AE(1,1)=1./ELEN 
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230 AE(1,2)=(-1. /ELEN) 

240 AE(2,1)=AE(1,2) 

250 AE(2,2)=AE(1,1) 

260 FS1EC1)=3. *ELEN 

270 FS1E(2)=FS1E(1) 

272 Fl=(ALPHA**4)*ELEN/2. 

274 F2=2. *(ALPHA**3)*(ELEN**2)/3. 

276 F3=(ALPHA**2)*(ELEN**3)/2. 

278 F4=ALPHA*(ELEN**4)/5. 

280 F5=(ELEN**5)/30. 

287 FS2E(1)=(-9. )*(F1 + F2 + F3 + F4 + F5) 

290 FS2E(2)=(-9)*(F1 + 2.*F2 + 3.*F3 + 4.*F4 + 5.*F5) 

300 DO 370 11=1,2 

310 DO 350 JJ=1,2 

320 IN=ICQRR(IEL,11) 

330 JN=ICORR(IEL,JJ) 

340 A(IN,JN)=A(IN,JN) - AE(II,JJ) 

350 CONTINUE 

360 FS(IN)=FS1E(II) + FS2E(II) + FS(IN) 

370 CONTINUE 

372 ALPHA=ALPHA + ELEN 

375 CONTINUE 

420 RETURN 

430 END 
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k 

k 


100 SUBROUTINE NU2KAI(IET) 

102 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
:UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),ITYPE,NEL,NSNP 

104 DIMENSION WKAREA(40600), DIF(IOO), BE(2,2), USTAR(IOO), FT(IOO) 

C SELECT METHOD OF DETERMINING USTAR 

105 PRINT*,’SELECT METHOD OF U* DETERMINATION. ' 

106 PRINT*,'1: U* = U' 

107 PRINT*,'2: U* = AVERAGE OF LAST TWO COMPUTED VALUES OF U' 

109 READC6,*) METHU 

C SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR 

116 PRINT*,'SELECT METHOD OF DETERMINING UNSTEADY B MATRIX. ' 

117 PRINT*,'1: MIDPOINT APPROXIMATION FOR U OVER THE ELEMENT. ’ 

119 PRINT*,'2: U LINEARIZED OVER THE LENGTH OF THE ELEMENT. ' 

120 READ(6,*) METHBM 

C CALL SUBROUTINE SETIME TO BEGIN TIMING ITERATION PROCESS 

121 CALL SETIME 

C BEGIN ITERATION PROCESS 

122 DO 450 ITER=1,200 

C DETERMINE VALUE OF U* AT EACH NODE 

130 DO 138 IU=1,NSNP 

131 IF (METHU.EQ. 1) THEN 

132 USTAR(IU)=U(IU) 

133 ELSE 

134 USTAR(IU)=(U(IU)+UOLD(IU))/2. 

137 ENDIF 

138 CONTINUE 

C DETERMINE UNSTEADY ELEMENT B MATRIX 

140 DO 210 IEL=1,NEL 

145 IF (METHBM. EQ. 1) THEN 

146 BE(1,1)=(ELEN/6. )*(USTAR(IEL)+USTAR(IEL+1)) 

147 BE(1,2)=(ELEN/12.)*(USTAR(IEL)+USTAR(IEL+1)) 

148 BE(2,1)=BE(1,2) 

149 BE(2,2)=BE(1,1) 

150 ELSE 

151 BE(1,1)=(ELEN/12. )*(3.*USTAR(IEL) + USTAR(IEL+1)) 

152 BE(1,2)=(ELEN/12. )*(USTAR(IEL) + USTAR(IEL+1)) 


C * SUBROUTINE NU2KAI 

C * 

C * THIS SUBROUTINE PERFORMS THE ITERATIVE SOLUTION PROCESS FOR 
C * MAIN PROGRAM NU2KA. 
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153 BE(2,1)=BE(1,2) 

154 BE(2,2)=(ELEN/12. )*(USTAR(IEL) + 3. *USTAR(IEL+1)) 

156 ENDIF 

170 DO 200 11=1,2 

175 DO 195 JJ=1,2 

180 IN=ICORR(IEL,II) 

185 JN=ICORR(IEL,JJ) 

190 B(IN,JN)=BE(II,JJ) + B(IN,JN) 

195 CONTINUE 

200 CONTINUE 

210 CONTINUE 

C DETERMINE TOTAL SYSTEM MATRIX 

220 DO 240 IP=1,NSNP 

221 DO 232 JP=1,NSNP 

230 C(IP,JP)=A(IP,JP)-B(IP,JP) 

C RESET B MATRIX TO ZERO AND LET UOLD=U AND FT=FS 

231 B(IP,JP)=0. 

232 CONTINUE 

235 UOLD(IP)=U(IP) 

236 FT(IP)=FS(IP) 

240 CONTINUE 

C IMPOSE BOUNDARY CONDITIONS 

241 C(l,l)=l. 

242 C(1,2)=0. 

243 FT(1)=ULBC 

244 IF (ITYPE.EQ. 1) THEN 

245 C(NSNP,NSNP-1)=0. 

246 C(NSNP,NSNP)=1. 

247 FT(NSNP)=URBC 

248 ELSE 

249 FT(NSNP)=FT(NSNP)-URBC 

250 ENDIF 

255 M=1 

260 IDGT=3 

270 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE SET OF LINEAR ALGEBRAIC EQUATIONS 

280 CALL LEQT2F(C,M,NSNP,IQ,FT,IDGT,WKAREA,IER) 

290 DO 310 NEW=1,NSNP 

300 U(NEW)=FT(NEW) 

C WRITE(*,*) 'UNEW=',U(NEW) 

C TEST FOR CONVERGENCE 

306 DIF(NEW)=ABS(U(NEW)-UOLD(NEW)) 

310 CONTINUE 

320 DIFMAX=DIF(1) 

325 NMAX=1 

330 DO 390 IJ=1,NEL 
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340 IF (DIF(IJ+1). GE.DIF(IJ)) THEN 

350 DIFMAX=DIF(IJ+1) 

355 NMAX=IJ+1 

360 ELSE 

370 CONTINUE 

380 ENDIF 

390 CONTINUE 

405 IF (ABS(DIFMAX/U(NMAX)). LT. CONV) THEN 

410 GO TO 460 

420 ELSE 

430 CONTINUE 

440 ENDIF 

450 CONTINUE 

460 CONTINUE 

C CALL SUBROUTINE GETIME TO OBTAIN CPU TIME FOR ITERATION PROCESS 


461 CALL GETIME(IET) 

C OUTPUT HEADER INFORMATION 


462 WRITEC6,464) 

463 WRITE(30,464) 

464 FORMATCIXEQUATION: U" - U**2 = 6 - 9X**4') 

465 IF (ITYPE.EQ. 1) THEN 

466 WRITE(6,468) COORD(1),ULBC,COORD(NSNP),URBC 

467 WRITEC30,468) COORDC1).ULBC.COORDCNSNP).URBC 

468 FORMATCIX,'B. C. : U(',F2. 0,')=',F2.0,*; U(',F2. 0,')=’,F4. 0,/) 

469 ELSE 

470 WRITE(6,472) COORD(1).ULBC,COORDCNSNP).URBC 

471 WRITE(30,472) COORDC1).ULBC,COORDCNSNP).URBC 

472 FORMATCIX,'B.C.: UC',F2. 0,')=',F2. 0,'; DU/DXC’,F2.0,’)=',F4.0,/) 

473 ENDIF 

475 IF CMETHU.EQ. 1) THEN 

476 WRITEC6.478) 

477 WRITEC30.478) 

478 FORMATCIX,'ITERATION METHOD: U*=U',/) 

479 ELSE 

480 WRITEC6.482) 

481 WRITEC 30,482) 

482 FORMATCIX,'ITERATION METHOD: U*=CU+UOLD)/2’,/) 

487 ENDIF 

488 IF CMETHBM.EQ. 1) THEN 

489 WRITEC6.491) 

490 WRITEC30,491) 

491 FORMATCIX,’METHOD OF B MATRIX INTEGRAL EVALUATION: MIDPOINT’,/) 

496 ELSE 

497 WRITEC6.499) 

498 WRITEC30,499) 

499 FORMATCIX,’METHOD OF B MATRIX INTEGRAL EVALUATION: LINEAR’,/) 

500 ENDIF 

502 IF CITER. GE. 200) THEN 

503 WRITEC6,505) 

504 WRITEC30,505) 

505 FORMATCIX,’CONVERGENCE NOT OBTAINED AFTER 200 ITERATIONS.’) 

506 ELSEIF CABSCUCNMAX)). GT. C10. **20). OR. ABSCUCNSNP-1)). GT. 
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:(10.**20)) THEN 

507 WRITE(6,509) 

508 WRITE(30,509) 

509 F0RMAT(IX,'SOLUTION PROCESS DIVERGES.') 

510 ELSE 

511 WRITE(6,520) ITER,NEL 

515 WRITE(30,520) ITER.NEL 

520 FORMAT(IX,'CONVERGENCE OBTAINED AFTER ',13,' ITERATIONS USING , 

:13,' ELEMENTS.',/) 

525 ENDIF 
530 RETURN 
540 END 








SUBROUTINE CLEXTA 


C 

c 

c 

c 

c 

c 


* 

* 


'if 


* 

* 


* THIS SUBROUTINE COMPUTES THE EXACT SOLUTION, U=3X**2, FOR * 

* MAIN PROGRAM NU2KA AT THE SPECIFIED NODAL POINTS. * 


100 SUBROUTINE CLEXTA 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
:UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),ITYPE,NEL,NSNP 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) = 3.*C00RD(NN)**2 

150 CONTINUE 
160 RETURN 
170 END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


* PROGRAM NU2KB * 

* THIS PROGRAM SOLVES THE NONLINEAR SECOND ORDER DIFFERENTIAL * 

* EQUATION: * 

* U" + U**2 = 60X + 100X**6; UEXACT=10X**3 WITH VARIABLE * 

* DOMAIN BY LINEARIZING THE U**2 TERM AS USTAR*U AND KEEPING IT * 

* ON THE LEFT SIDE OF THE EQUATION. THE USER SELECTS: * 

* 1) NUMBER OF ELEMENTS * 

* 2) SIZE OF DOMAIN * 

* 3) X AND U(X) AT THE LEFT BOUNDARY * 

* 4) U(X) OR U'(X) AT THE RIGHT BOUNDARY * 

* 5) ITERATION STRATEGY FOR DETERMINING U* * 

* 6) APPROXIMATION TECHNIQUE FOR THE EXCITATION INTEGRAL * 


110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),ITYPE,NEL,NSNP 
115 CONV=. 0001 


C READ IN PARAMETERS FROM DATA FILE 

130 READ(29,*) NEL,TLEN,COORD(1),ULBC,ITYPE,URBC 

C CALCULATE NUMBER OF NODAL POINTS 


135 NSNP=NEL+1 


C DETERMINE ELEMENT SIZE OF EQUAL LENGTHS 
137 ELEN=TLEN/FLOAT(NEL) 

C ESTABLISH LOCAL TO GLOBAL CORRESPONDENCE AND X COORDINATE OF 
C EACH NODE 

160 DO 169 IEL=1,NEL 

162 ICORR(IEL,1)=IEL 

163 IC0RR(IEL,2)=IEL+1 

164 C00RD(IEL+1)=C00RD(IEL)+ELEN 

169 CONTINUE 

C CALL SUBROUTINE NU2KBM TO CREATE A MATRIX AND F VECTOR 

170 CALL NU2KBM 

C CALL SUBROUTINE NU2KBI TO PERFORM SOLUTION ITERATION 

180 CALL NU2KBI(IET) 

C CALL SUBROUTINE CLEXTB TO COMPUTE EXACT SOLUTION U=3X**2 

190 CALL CLEXTB 

C CALL SUBROUTINE CLOTPT TO PRINT OUT DATA, COMPUTATIONAL EFFICIENCY 
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200 CALL CLOTPT(CPUSTAR,IET) 
210 END 






c 

c 

c 

c 

c 

c 


SUBROUTINE NU2KBM 


* 

* 

* THIS SUBROUTINE COMPUTES THE A MATRIX AND F VECTOR FOR MAIN 

* PROGRAM NU2KB. 


* 

* 

* 

* 


100 SUBROUTINE NU2KBM 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,IC0RR(100,2),ITYPE,NEL,NSNP 
120 DIMENSION AE(2,2), FS1E(2), FS2E(2) 

C ZERO OUT A MATRIX AND ALL VECTORS 


140 DO 210 IZ = 1,NSNP 
150 FS(IZ) = 0. 

154 U(IZ) =SQRT(60. *COORD(IZ) + 100.*COORD(IZ)**6) 

156 UOLD(IZ)=0. 

160 DO 200 JZ = 1,NSNP 

170 A(IZ,JZ) = 0. 

175 B(IZ,JZ) = 0. 

176 C(IZ,JZ) = 0. 

200 CONTINUE 

210 CONTINUE 

C ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 


213 ALPHA=0. 

215 DO 375 IEL=1,NEL 
220 AE(1, 1)=1. /ELEN 

230 AE(1,2)=(-1. /ELEN) 

240 AE(2,1)=AE(1,2) 

250 AE(2,2)=AE(1,1) 

260 FS1E(1)=30.*ALPHA*ELEN + 10.*ELEN**2 

270 FS1E(2)=30. *ALPHA*ELEN + 20. *ELEN**2 

272 F1=50.*(ALPHA**6)*ELEN 

2 74 F2=100. *(ALPHA**5)*(ELEM**2) 

276 F3=125.*(ALPHA**4)*(ELEN**3) 

278 F4=100.*(ALPHA**3)*(ELEN**4) 

280 F5=50.*(ALPHA**2)*(ELEN**5) 

281 F6-100. *ALPHA*(ELEN**6)/7. 

282 F7=25.*(ELEN**7)/16. 

287 FS2E(1)=F1 + F2 + F3 + F4 + F5 + F6 + F7 

290 FS2E(2)=F1 + 2.*F2 + 3.*F3 + 4.*F4 + 5.*F5 + 6.*F6 + 7*F7 

300 DO 370 11=1,2 

310 DO 350 JJ=1,2 

320 IN=ICORR(IEL,II) 

330 JN=ICORR(IEL,JJ) 

340 A(IN,JN)=A(IN,JN) - AE(II,JJ) 

350 CONTINUE 

360 FS(IN)=FS1E(II) + FS2E(II) + FS(IN) 

370 CONTINUE 

372 ALPHA=ALPHA + ELEN 

375 CONTINUE 
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420 RETURN 
430 END 
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c 

c 

c 

c 

c 

c 


SUBROUTINE NU2KBI 


* 

★ 

* THIS SUBROUTINE PERFORMS THE ITERATIVE SOLUTION PROCESS FOR 

* MAIN PROGRAM NU2KA. 


* 

★ 

* 

* 


100 SUBROUTINE NU2KBI(IET) 

102 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 

: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

:TLEN,ICORRC100,2),ITYPE.NEL.NSNP 

104 DIMENSION WKAREA(40600), DIF(IOO), BE(2,2), USTAR(IOO), FT(IOO) 

C SELECT METHOD OF DETERMINING USTAR 

105 PRINT*,'SELECT METHOD OF U* DETERMINATION.' 

106 PRINT*,'1: U* = U' 

107 PRINT*,'2: U* = AVERAGE OF LAST TWO COMPUTED VALUES OF U' 

108 PRINT*,'3: U* = WEIGHTED AVERAGE OF LAST TWO COMPUTED VALUES OF U' 

109 READ(6,*) METHU 

110 IF (METHU. EQ. 3) THEN 

111 PRINT*,'CHOOSE WEIGHTING VALUES A AND B' 

112 READ(6,*) AW,BW 

113 ELSE 

114 CONTINUE 

115 ENDIF 


C SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR 

116 PRINT*,'SELECT METHOD OF DETERMINING UNSTEADY B MATRIX. ’ 

117 PRINT*,'1: MIDPOINT APPROXIMATION FOR U OVER THE ELEMENT. ' 

119 PRINT*,'2: U LINEARIZED OVER THE LENGTH OF THE ELEMENT.’ 

120 READC6,*) METHBM 

C CALL SUBROUTINE SETIME TO BEGIN TIMING ITERATION PROCESS 


121 CALL SETIME 


C BEGIN ITERATION PROCESS 


122 DO 450 ITER=1,200 
C DETERMINE VALUE OF U* AT EACH NODE 


130 DO 138 IU=1,NSNP 

131 IF (METHU.EQ.1) THEN 

132 USTARCIU)=U(IU) 

133 ELSEIF (METHU. EQ. 2) THEN 

134 USTARCIU)=(U(IU)+U0LD(IU))/2. 

135 ELSE 

136 USTARCIU)=(AW*U(IU)+BW*UOLD(IU))/(AW+BW) 

137 ENDIF 

138 CONTINUE 

C DETERMINE UNSTEADY ELEMENT B MATRIX 
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140 DO 210 IEL=1,NEL 

145 IF (METHBM.EQ. 1) THEN 

146 BE(1,1)=(ELEN/6. )*(USTAR(IEL)+USTAR(IEL+1)) 

147 BE(1,2)=(ELEN/12. )*(USTAR(IEL)+USTAR(IEL+1)) 

148 BE(2,1)=BE(1,2) 

149 BE(2,2)=BE(1,1) 

150 ELSE 

151 BE( 1,1)=(ELEN/12. )*(3. *USTAR(IEL) + USTAR(IEL+1)) 

152 BE(1,2)=(ELEN/12. )*(USTAR(IEL) + USTAR(IEL+1)) 

153 BE(2,1)=BE(1,2) 

154 BE(2,2)=(ELEN/12. )*(USTAR(IEL) + 3. *USTAR(IEL+1)) 

156 ENDIF 

170 DO 200 11=1,2 

175 DO 195 JJ=1,2 

180 IN=ICORR(IEL,II) 

185 JN=ICORR(IEL,JJ) 

190 B(IN,JN)=BE(II,JJ) + B(IN.JN) 

195 CONTINUE 

200 CONTINUE 

210 CONTINUE 


C DETERMINE TOTAL SYSTEM MATRIX 

220 DO 240 IP=1,NSNP 

221 DO 232 JP=1,NSNP 

230 C(IP,JP)=A(IP,JP)+B(IP,JP) 

C RESET B MATRIX TO ZERO AND LET UOLD=U AND FT=FS 

231 B(IP,JP)=0. 

232 CONTINUE 

235 UOLD(IP)=U(IP) 

236 FT(IP)=FS(IP) 

240 CONTINUE 


C IMPOSE BOUNDARY CONDITIONS 

241 C(l,l)=l. 

242 C(l,2)=0. 

243 FT(1)=ULBC 

244 IF (ITYPE.EQ. 1) THEN 

245 C(NSNP,NSNP-1)=0. 

246 C(NSNP,NSNP)=1. 

247 FT(NSNP)=URBC 

248 ELSE 

249 FT(NSNP)=FT(NSNP)-URBC 

250 ENDIF 

255 M=1 

260 IDGT=3 

270 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE SET OF LINEAR ALGEBRAIC EQUATIONS 

280 CALL LEQT2F(C,M,NSNP,IQ,FT,IDGT,WKAREA,IER) 

290 DO 310 NEW=1,NSNP 
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300 U(NEW)=FT(NEW) 

C WRITEC*,*) 'UNEW=',U(NEW) , 

C TEST FOR CONVERGENCE 

306 DIF(NEW)=ABS(U(NEW)"UOLD(NEW)) 

310 CONTINUE 

320 DIFMAX=DIF(1) 

325 NMAX=1 

330 DO 390 IJ=1,NEL 

340 IF (DIF(IJ+1). GE. DIF(IJ)) THEN 

350 DIFMAX=DIF(IJ+1) 

355 NMAX=IJ+1 

360 ELSE 

370 CONTINUE 

380 ENDIF 

390 CONTINUE 

405 IF (ABS(DIFMAX/U(NMAX)). LT. CONV) THEN 

410 GO TO 460 

420 ELSE 

430 CONTINUE 

440 ENDIF 

450 CONTINUE 

460 CONTINUE 

C CALL SUBROUTINE GETIME TO OBTAIN CPU TIME FOR ITERATION PROCESS 

461 CALL GETIME(IET) 

C OUTPUT HEADER INFORMATION 

462 WRITE(6,464) 

463 WRITEC30,464) 

464 FORMATCIX,'EQUATION: U" + U**2 = 60X + 100X**6') 

465 IF (ITYPE.EQ. 1) THEN 

466 WRITE(6,468) COORD(1),ULBC,COORD(NSNP),URBC 

467 WRITE(30,468) COORD(1),ULBC.COORD(NSNP),URBC 

468 FORMAT(IX,'B.C.: U(',F2. 0,')=’,F3.0,'; U(',F2. 0,')=’,F5.0,/) 

469 ELSE 

470 WRITE(6,472) COORD(1),ULBC,COORD(NSNP),URBC 

471 WRITE(30,472) COORD(1),ULBC,COORD(NSNP),URBC 

472 FORMATCIX,'B.C.: UC',F2. 0,',F3. 0,';DU/DXC',F2.0,')«',F5.0,/) 

473 ENDIF 

475 IF CMETHU.EQ. 1) THEN 

476 WRITEC6.478) 

477 WRITEC 30,478) 

478 FORMATCIX,'ITERATION METHOD: U*=U', / ) 

479 ELSEIF CMETHU.EQ. 2) THEN 

480 WRITEC6.482) 

481 WRITEC30,482) 

482 FORMATCIX,'ITERATION METHOD: U*=(U+UOLD)/2 1 ,/) 

483 ELSE 

484 WRITEC6,486)AW,BW,AW,BW 

485 WRITEC30,486)AW,BW,AW,BW 

486 FORMATCIX,'ITERATION METHOD: U*=(',F3.0,'*U +',F3.0,'*UOLD)/C’, 
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:F3.0,'+',F3.0,')',/) 

487 ENDIF 

488 IF (METHBM.EQ. 1) THEN 

489 WRITE(6,491) 

490 WRITE(30,491) 

491 F0RMAT(IX,'METHOD OF B MATRIX INTEGRAL EVALUATION: MIDPOINT',/) 

496 ELSE 

497 WRITE(6,499) 

498 WRITE(30,499) 

499 FORMAT(IX,'METHOD OF B MATRIX INTEGRAL EVALUATION: LINEAR’,/) 

500 ENDIF 

502 IF (ITER. GE. 200) THEN 

503 WRITE(6,505) 

504 WRITE(30,505) 

505 FORMAT(IX,'CONVERGENCE NOT OBTAINED AFTER 200 ITERATIONS.’) 

506 ELSEIF (ABS(U(NMAX)). GT. (10. **20). OR. ABS(U(NSNP-1)). GT. 

:(10.**20)) THEN 

507 WRITE(6,509) 

508 WRITE(30,509) 

509 FORMAT(IX,'SOLUTION PROCESS DIVERGES. ') 

510 ELSE 

511 WRITE(6,520) ITER.NEL 

515 WRITE(30,520) ITER.NEL 

520 FORMAT(IX,'CONVERGENCE OBTAINED AFTER ',13,' ITERATIONS USING ', 

:13,' ELEMENTS. ',/) 

525 ENDIF 
530 RETURN 
540 END 
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c 

c 

c 

c 

c 

c 


* AA A A fr AA 'A’An fr **"* * A AAcAAA ' A ;Ar*ifcA;Ari V AAA kiviti r hi 'f / r ifi t iri r /tiriri r / r / r k / t i ri r / r iti r ki r f r / f kAI r k-h'kr/t 

* SUBROUTINE CLEXTB * 

* * 


* THIS SUBROUTINE COMPUTES THE EXACT SOLUTION, U=10X**3, FOR * 

* MAIN PROGRAM NU2KB AT THE SPECIFIED NODAL POINTS. * 


100 SUBROUTINE CLEXTB 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORRC100,2),ITYPE,NEL,NSNP 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) = 10. *COORD(NN)**3 

150 CONTINUE 
160 RETURN 
170 END 
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Q it * *** * * ** * ******** ** ** **. *** A A* ** * *********iWf*** **AA» *********^*** 

C * SUBROUTINE CLOTPT * 

C * THIS SUBROUTINE COMPUTES THE PER CENT ERROR BETWEEN THE EXACT * 

C * AND FEM SOLUTIONS, CPU* FOR THE ITERATION PROCESS, AND PRINTS * 

C * OUT ALL DATA IN TABULAR FORM FOR MAIN PROGRAMS NU2KA & NU2KB^* 

100 SUBROUTINE CLOTPT(CPUSTAR.IET) 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 

: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: TLEN,ICORR(100,2),ITYPE,NEL,NSNP 
115 SUMDIF=0. 

C CALCULATE PER CENT ERROR AT EACH NODE AND SUM THE ABSOLUTE VALUE 
C OF ALL THE ERRORS 

120 DO 150 IK=2,NSNP 

130 UDIF(IK)=100. *(U(IK)-UEXT(IK))/UEXT(IK) 

140 SUMDIF=SUMDIF + ABS(UDIF(IK)) 

150 CONTINUE 

160 UDIF(1)=U(1)-UEXT(1) 

C COMPUTE THE ELAPSED TIME OF THE ITERATION PROCESS 

164 ELTIME=IET*. 000026 

165 WRITE(6,169) ELTIME 

166 WRITE(30,169)ELTIME 

169 FORMAT(IX,'ELAPSED TIME FOR THE ITERATION PROCESS IS ',F9.4, 

: ' SECONDS.') 

C OUTPUT DATA IN TABULAR FORMAT 

170 WRITE(6,180) 

175 WRITE(30,180) 

180 FORMAT(/,IX,'X-COORD',3X,'U EXACT',3X,’U FEM',7X,'% DIFF') 

190 WRITE(6,200) (COORD(NP),UEXT(NP),U(NP),UDIF(NP), NP=1,NSNP) 

195 WRITE(30,200) (COORD(NP),UEXT(NP),U(NP),UDIF(NP), NP=1,NSNP) 

200 FORMAT(/,2X,F5.3,4X,F9.4,3X,F9. 4,4X,F5. 1) 

C CALCULATE CPU* FOR THE ITERATION PROCESS 

205 CPUSTAR=ELTIME*SUMDIF/NSNP 
210 WRITE(6,220) CPUSTAR 

215 WRITE(30,220) CPUSTAR 

220 F0RMAT(/,1X, CPU* FOR THE ITERATION PROCESS IS ',F9.4,' SECONDS.') 
230 RETURN 

240 END 
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APPENDIX G. PROGRAM LISTINGS FOR QUASILINEARIZATION 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


* PROGRAM NU2QA * 

* THIS PROGRAM SOLVES THE NONLINEAR SECOND ORDER DIFFERENTIAL * 

* EQUATION: * 

* U" - U**2 = 6 - 9X**4; UEXACT=3X**2 WITH VARIABLE DOMAIN * 

* BY THE PROCESS OF QUASILINEARIZATION. THE USER SELECTS: * 

* 1) NUMBER OF ELEMENTS * 

* 2) SIZE OF DOMAIN * 

* 3) X AND U(X) AT THE LEFT BOUNDARY * 

* 4) U(X) OR U'(X) AT THE RIGHT BOUNDARY * 

* 5) ITERATION STRATEGY FOR DETERMINING U* * 

* 6) INTERPOLATION STRATEGY FOR THE B MATRIX INTEGRAL * 

* 7) INTERPOLATION STRATEGY FOR THE EXCITATION INTEGRAL * 


110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: ICORR(100,2),ITYPE,NEL,NSNP 
115 CONV=. 0001 


C READ IN PARAMETERS FROM DATA FILE 

130 READ(29,*) NEL.TLEN,COORD(1),ULBC,ITYPE.URBC 

C CALCULATE NUMBER OF NODAL POINTS 


135 NSNP=NEL+1 

C DETERMINE ELEMENT SIZE OF EQUAL LENGTHS 

137 ELEN=TLEN/FLOAT(NEL) 

C ESTABLISH LOCAL TO GLOBAL CORRESPONDENCE AND X COORDINATE OF 
C EACH NODE 


160 DO 169 IEL=1,NEL 

162 IC0RR(IEL,1)=IEL 

163 IC0RR(IEL,2)=IEL+1 

164 C00RD(IEL+1)=C00RD(IEL)+ELEN 

169 CONTINUE 

C CALL SUBROUTINE NU2QAM TO CREATE A MATRIX AND F VECTOR 

170 CALL NU2QAM 

C CALL SUBROUTINE NU2QAI TO PERFORM SOLUTION ITERATION 

180 CALL NU2QAI(IET) 

C CALL SUBROUTINE QLEXTA TO COMPUTE EXACT SOLUTION U=3X**2 
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190 CALL QLEXTA 

C CALL SUBROUTINE QLOTPT TO PRINT OUT DATA, COMPUTATIONAL EFFICIENCY 

200 CALL QLOTPT(CPUSTAR,IET) 

210 END 





c 

c 

c 

c 

c 

c 



* THIS SUBROUTINE COMPUTES THE A MATRIX AND STEADY F VECTOR FOR * 

* MAIN PROGRAM NU2QA. * 


100 SUBROUTINE NU2QAM 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(ICO), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

:ICORR(100,2),ITYPE,NEL,NSNP 
120 DIMENSION AE(2,2), FS1E(2), FS2E(2) 

C ZERO OUT A MATRIX AND ALL VECTORS 


140 DO 210 IZ = l.NSNP 
150 FS(IZ) = 0. 

154 U(IZ) =SQRT(ABS(9.*C00RD(IZ)**4-6.)) 

156 UOLD(IZ)=0. 

160 DO 200 JZ = l.NSNP 

170 A(IZ,JZ) = 0. 

175 B(IZ,JZ) = 0. 

176 C(IZ,JZ) = 0. 

200 CONTINUE 

210 CONTINUE 


C ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 


213 

ALPHA=0. 


215 

DO 375 IEL=1,NEL 


220 

AE(1,1)=1./ELEN 


230 

AE(1,2)=(-1. /ELEN) 


240 

AE(2,1)=AE(1,2) 


250 

AE(2,2)=AE(1,1) 


260 

FS1E(1)=3. *ELEN 


270 

FS1E(2)=FS1E(1) 


272 

Fl=(ALPHA**4)*ELEN/2. 


274 

F2=2. *(ALPHA**3)*(ELEN**2)/3. 


276 

F3=(ALPHA**2)*(ELEN**3)/2. 


278 

F4=ALPHA*(ELEN**4)/5. 


280 

F5=(ELEN**5)/30. 


287 

FS2E(l)=(-9. )*(F1 + F2 + F3 + F4 

+ F5) 

290 

FS2E(2)=(-9)*(F1 + 2.*F2 + 3.*F3 

+ 4.*F4 

300 

DO 370 11=1,2 


310 

DO 350 JJ=1,2 


320 

IN=ICORR(IEL.II) 


330 

JN=ICORR(IEL,JJ) 


340 

A(IN,JN)=A(IN,JN) - AE(II,JJ) 

350 

CONTINUE 


360 

FS(IN)=FS1E(II) + FS2E(II) + FS(IN) 

370 

CONTINUE 


372 

ALPHA=ALPHA + ELEN 


375 

CONTINUE 


420 

RETURN 


430 

END 





SUBROUTINE NU2QAI 


C 

c 

c 

c 

c 

c 


■kirk********** A * ****** 
* 

* 



* THIS SUBROUTINE PERFORMS THE ITERATIVE SOLUTION PROCESS FOR 

* MAIN PROGRAM NU2QA. 


100 SUBROUTINE NU2QAI(IET) 

102 COMMON A( 100,100) ,FS( 100) ,B( 100,100) ,C( 100, lOO'l .U( 100) ,UOLD( 100) , 

: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: ICORR(100,2),ITYPE,NEL,NSNP 

104 DIMENSION WKAREA(40600), DIF(IOO), BE(2,2), USTAR(IOO), FT(IOO), 

:FUE(2),FU(100) 

C SELECT METHOD OF DETERMINING USTAR 

105 PRINT*,’SELECT METHOD OF U* DETERMINATION. ’ 

106 PRINT*,’1: U* = U* 

107 PRINT*,’2: U* = AVERAGE OF LAST TWO COMPUTED VALUES OF U* 

108 PRINT*,’3: U* = WEIGHTED AVERAGE OF LAST TWO COMPUTED VALUES OF U’ 

109 READ(6,*) METHU 

110 IF (METHU. EQ. 3) THEN 

111 PRINT*,’CHOOSE WEIGHTING VALUES A AND B’ 

112 READ(6,*) AW,BW 

113 ELSE 

114 CONTINUE 

115 ENDIF 


C SELECT METHOD OF DETERMINING UNSTEADY B MATRIX 

116 PRINT*,’SELECT METHOD OF DETERMINING UNSTEADY B MATRIX. ’ 

117 PRINT*,’1: MIDPOINT APPROXIMATION FOR U OVER THE ELEMENT. ’ 

119 PRINT*,’2: U LINEARIZED OVER THE LENGTH OF THE ELEMENT. ’ 

120 READ(6,*) METHBM 

C SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR 

121 PRINT*,’SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR. ’ 

122 PRINT*,’1: MIDPOINT APPROXIMATION’ 

123 PRINT*,’2: 1/4 - 3/4 APPROXIMATION’ 

124 PRINT*,’3: LINEAR’ 

125 READ(6,*) METHFU 

C CALL SUBROUTINE SETIME TO BEGIN TIMING ITERATION PROCESS 


127 CALL SETIME 


C BEGIN ITERATION PROCESS 


128 DO 450 ITER=1,200 

C DETERMINE VALUE OF U* AT EACH NODE AND SET VALUE OF UNSTEADY 
C FORCE VECTOR TO ZERO 


129 DO 138 IU=1,NSNP 

130 FU(IU)=0. 

131 IF (METHU.EQ. 1) THEN 
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132 USTAR(IU)=U(IU) 

133 ELSEIF (METHU.EQ. 2) THEN 

134 USTAR(IU)=(U(IU)+UOLD(IU))/2. 

135 ELSE 

136 USTAR(IU)=(AW*U(IU)+BW*UOLD(IU))/(AW+BW) 

137 ENDIF 

138 CONTINUE 

C DETERMINE UNSTEADY ELEMENT B MATRIX 

140 DO 210 IEL=1,NEL 

145 IF (METHBM.EQ. 1) THEN 

146 BE(l,l)=(ELEN/3. )*(USTAR(IEL)+USTAR(IEL+1)) 

147 BE(1,2)=(ELEN/6. )*(USTAR(IEL)+USTAR(IEL+1)) 

148 BE(2,1)=BE(1,2) 

149 BE(2,2)=BE(1,1) 

150 ELSE 

151 BE(1,1)=(ELEN/6. )*(3. *USTAR(IEL) + USTAR(IEL+1)) 

152 BE(l,2)=(ELEN/6. )*(USTAR(IEL) + USTAR(IEL+1)) 

153 BE(2,1)=BE(1,2) 

154 BE(2,2)=(ELEN/6. )*(USTAR(IEL) + 3.*USTAR(IEL+1)) 

156 ENDIF 

C DETERMINE SYSTEM B MATRIX BY DISTRIBUTING ELEMENT B MATRICES 

C ACCORDING TO THE LOCAL TO GLOBAL CORRESPONDENCE 

157 DO 163 11=1,2 

158 DO 162 JJ=1,2 

159 IN=ICORR(IEL,II) 

160 JN=ICORR(IEL,JJ) 

161 B(IN,JN)=BE(II,JJ) + B(IN,JN) 

162 CONTINUE 

163 CONTINUE 

C DETERMINE UNSTEADY ELEMENT FORCE VECTOR 

165 IF (METHFU.EQ. 1) THEN 

166 FUE(l)=(ELEN/2. )*((USTAR(IEL)+USTAR(IEL+1))/2. )**2 

167 FUE(2)=FUE(1) 

169 ELSEIF (METHFU.EQ. 2) THEN 

171 FUE(l)=(ELEN/2. )*(3.*USTAR(IEL)/4. + USTAR(IEL+1)/4. )**2 

172 FUE(2)=(ELEN/2. )*(USTAR(IEL)/4. + 3.*USTAR(IEL+1)/4. )**2 

173 ELSE 

174 FUE(l)=ELEN*(USTAR(IEL)**2/4. +USTAR(IEL)*USTAR(IEL+l)/6. 

: + USTAR(IEL+1)**2/12. ) 

175 FUE(2)=ELEN*(USTAR(IEL)**2/12. + USTAR(IEL)*USTAR(IEL+1)/6. 

+ USTAR(IEL+l)**2/4.) 

176 ENDIF 

C DETERMINE UNSTEADY SYSTEM FORCE VECTOR BY DISTRIBUTING ELEMENTAL 

C FORCE VECTORS ACCORDING TO THE LOCAL TO GLOBAL CORRESPONDENCE 

177 DO 180 11=1,2 

178 IN=ICORR(IEL,II) 

179 FU(IN)=FUE(II) + FU(IN) 
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180 CONTINUE 

210 CONTINUE 

C DETERMINE TOTAL SYSTEM MATRIX 

220 DO 240 IP=1,NSNP 

221 DO 232 JP=1,NSNP 

230 C(IP,JP)=A(IP,JP)-B(IP,JP) 

C RESET B MATRIX TO ZERO 

231 B(IP,JP)=0. 

232 CONTINUE 

C UPDATE VALUE OF U AT THE PREVIOUS ITERATION 

235 UOLD(IP)=U(IP) 

C DETERMINE TOTAL SYSTEM FORCE VECTOR 

236 FT(IP)=FS(IP)-FU(IP) 

240 CONTINUE 

C IMPOSE BOUNDARY CONDITIONS 

241 C(l,l)=l. 

242 C(l,2)=0. 

243 FT(1)=ULBC 

244 IF (ITYPE.EQ. 1) THEN 

245 C(NSNP,NSNP-1)=0. 

246 C(NSNP,NSNP)=1. 

247 FT(NSNP)=URBC 

248 ELSE 

249 FT(NSNP)=FT(NSNP)-URBC 

250 ENDIF 

255 M=1 

260 IDGT=3 

270 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE SET OF LINEAR ALGEBRAIC EQUATIONS 

280 CALL LEQT2F(C,M,NSNP,IQ,FT,IDGT,WKAREA,IER) 

290 DO 310 NEW=1,NSNP 

300 U(NEW)=FT(NEW) 

C WRITE(*,*) 1 UNE¥=',U(NEW) 

C TEST FOR CONVERGENCE 

306 DIF(NEW)=ABS(U(NEW)-UOLD(NEW)) 

310 CONTINUE 

320 DIFMAX=DIF(1) 

325 NMAX=1 

330 DO 390 IJ=1,NEL 

340 IF (DIF(IJ+1). GE. DIF(IJ)) THEN 

350 DIFMAX=DIF(IJ+1) 

355 NMAX=IJ+1 
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360 ELSE 

370 CONTINUE 

380 ENDIF 

390 CONTINUE 

400 IF (U(NMAX).EQ.O. ) GO TO 450 

405 IF (ABS(DIFMAX/U(NMAX)). LT. CONV) THEN 

410 GO TO 460 

420 ELSE 

430 CONTINUE 

440 ENDIF 

450 CONTINUE 

460 CONTINUE 

C CALL SUBROUTINE GETIME TO OBTAIN CPU TIME FOR ITERATION PROCESS 

461 CALL GETIME(IET) 

C OUTPUT HEADER INFORMATION 

462 WRITE(6,464) 

463 WRITE(30,464) 

464 FORMATCIX,'EQUATION: U" - U**2 = 6 - 9X**4') 

465 IF (ITYPE.EQ. 1) THEN 

466 WRITE(6,468) COORD(1),ULBC,COORD(NSNP),URBC 

467 WRITEC30,468) COORD(1),ULBC.COORD(NSNP),URBC 

468 F0RMAT(1X,'B.C. : U(',F2. 0,')=',F2. 0,'; U(',F2. 0,')=',F4. 0,/) 

469 ELSE 

470 WRITE(6,472) COORD(1),ULBC,COORD(NSNP),URBC 

471 WRITE(30,472) COORD(1),ULBC.COORD(NSNP),URBC 

472 FORMATCIX, ’ B. C. : U( 1 ,F2.0,')=’,F2.0,’;DU/DX(’,F2.0,’)=’,F4. 0,/) 

473 ENDIF 

475 IF (METHU.EQ. 1) THEN 

476 WRITE(6,478) 

477 WRITE(30,478) 

478 FORMATCIX,'ITERATION METHOD: U*=U',/) 

479 ELSEIF (METHU.EQ. 2) THEN 

480 WRITE(6,482) 

481 WRITE(30,482) 

482 FORMATCIX,'ITERATION METHOD: U*=(U+UOLD)/2',/) 

483 ELSE 

484 WRITE(6,486)AW,BW,AW,BW 

485 WRITEC30,486)AW,BW,AW,BW 

486 FORMAT(IX,'ITERATION METHOD: U*=(',F4. 1,'*U +',F4.1,'*UOLD)/(', 

:F4.1,'+',F4.1, ’)',/) 

487 ENDIF 

488 IF (METHBM.EQ. 1) THEN 

489 WRITE(6,491) 

490 WRITEC30,491) 

491 FORMATCIX,'METHOD OF B MATRIX INTEGRAL EVALUATION: MIDPOINT',/) 

496 ELSE 

497 WRITE(6,499) 

498 WRITEC30,499) 

499 FORMATCIX,'METHOD OF B MATRIX INTEGRAL EVALUATION: LINEAR',/) 

500 ENDIF 

501 IF (METKFU.EQ. 1) THEN 

502 WRITEC6,504) 
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503 WRITE(30,504) 

504 F0RMAT(IX,'METHOD OF EXCITATION INTEGRAL EVALUATION: MIDPOINT',/) 

505 ELSEIF (METHFU. EQ. 2) THEN 

506 WRITE(6,508) 

507 WRITE(30,508) 

508 FORMAT(IX,'METHOD OF EXCITATION INTEGRAL EVALUATION: 1/4-3/4’,/) 

509 ELSE 

510 WRITE(6,512) 

511 WRITE(30,512) 

512 FORMATCIX,'METHOD OF EXCITATION INTEGRAL EVALUATION: LINEAR’,/) 

513 ENDIF 

515 IF (ITER. GE. 200) THEN 

516 WRITE(6,518) 

517 WRITE(30,518) 

518 FORMAT(IX,'CONVERGENCE NOT OBTAINED AFTER 200 ITERATIONS. ') 

519 ELSEIF (ABS(U(NMAX)).GT. (10. **20). OR. ABS(U(NSNP-1)). GT. 

:(10.**20)) THEN 

520 WRITE(6,522) 

521 WRITE(30,522) 

522 FORMAT(IX,'SOLUTION PROCESS DIVERGES.') 

523 ELSE 

524 WRITE(6,526) ITER.NEL 

525 WRITE(30,526) ITER.NEL 

526 FORMATCIX,'CONVERGENCE OBTAINED AFTER ',13,' ITERATIONS USING ', 
: 13,' ELEMENTS. ',/) 

530 ENDIF 
540 RETURN 
550 END 
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c 

c 

c 

c 

c 

c 


*** * * ** * * * * ****** * **************** A* * *A A* A****** * ** *** * ** * ******** 

* SUBROUTINE QLEXTA * 

* * 


* THIS SUBROUTINE COMPUTES THE EXACT SOLUTION, U=3X**2, FOR 

* MAIN PROGRAM NU2QA AT THE SPECIFIED NODAL POINTS. 


100 SUBROUTINE QLEXTA 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
:UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: ICORR(100,2),ITYPE,NEL,NSNP 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) = 3. *C00RD(NN)**2 

150 CONTINUE 
160 RETURN 
170 END 
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•) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


* PROGRAM NU2QB * 

* THIS PROGRAM SOLVES THE NONLINEAR SECOND ORDER DIFFERENTIAL * 

* EQUATION: * 

* u" + U**2 = 6OX + 100X**6; UEXACT=10X**3 WITH VARIABLE * 

* DOMAIN BY THE PROCESS OF QUASILINEARIZATION. THE USER SELECTS:* 

* 1) NUMBER OF ELEMENTS * 

* 2) SIZE OF DOMAIN * 

* 3) X AND U(X) AT THE LEFT BOUNDARY * 

* 4) U(X) OR U'(X) AT THE RIGHT BOUNDARY * 

* 5) ITERATION STRATEGY FOR DETERMINING U* * 

* 6) INTERPOLATION STRATEGY FOR THE B MATRIX INTEGRAL * 

* 7) INTERPOLATION STRATEGY FOR THE EXCITATION INTEGRAL * 


110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
:UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: ICORR(100,2),ITYPE,NEL,NSNP 
115 CONV=. 0000001 

C READ IN PARAMETERS FROM DATA FILE 
130 READ(29,*) NEL,TLEN,COORD(1),ULBC,ITYPE,URBC 

C CALCULATE NUMBER OF NODAL POINTS 
135 NSNP=NEL+1 

C DETERMINE ELEMENT SIZE OF EQUAL LENGTHS 
137 ELEN=TLEN/FLOAT(NEL) 

C ESTABLISH LOCAL TO GLOBAL CORRESPONDENCE AND X COORDINATE OF 
C EACH NODE 

160 DO 169 IEL=1,NEL 

162 ICORR(IEL,1)=IEL 

163 ICORR(IEL,2)=IEL+l 

164 COORD(IEL+l)=COORD(IEL)+ELEN 

169 CONTINUE 

C CALL SUBROUTINE NU2QBM TO CREATE A MATRIX AND F VECTOR 

170 CALL NU2QBM 

C CALL SUBROUTINE NU2QBI TO PERFORM SOLUTION ITERATION 

180 CALL NU2QBI(IET) 

C CALL SUBROUTINE QLEXTB TO COMPUTE EXACT SOLUTION U=3X**2 

190 CALL QLEXTB 

C CALL SUBROUTINE QLOTPT TO PRINT OUT DATA, COMPUTATIONAL EFFICIENCY 
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200 CALL QLOTPTC CPUSTAR,IET) 
210 END 







Q *- ki t*'kirk- ir k'kirkirk*i rk -kirk** rk ' k **irk-k-k**'k-kkkkkk A 'A' * ** * **** ** **** * * **** *& '*■'* 

C * SUBROUTINE NU2QBM * 

C * * 

C * THIS SUBROUTINE COMPUTES THE A MATRIX AND STEADY F VECTOR FOR * 

C * MAIN PROGRAM NU2QB. * a,uiaaaaaaIa 

100 SUBROUTINE NU2QBM 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
: UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: ICORR(100,2),ITYPE,NEL,NSNP 

120 DIMENSION AE(2,2), FS1E(2), FS2E(2), Gl(100), G2(100), G3(100) 

C ZERO OUT A MATRIX AND ALL VECTORS 

140 DO 210 IZ = l.NSNP 
145 FS(IZ) = 0. 

C IF (COORD(IZ).LE. 1. ) THEN 

C U(IZ)=0. 

C ELSE 

C Gl(IZ) = 60.*COORD(IZ) + 100.*COORD(IZ)**6 

C G2(IZ) = 1500. *COORD(IZ)**4/SQRT(G1(IZ)) 

C G3(IZ) = (60. +600. *COORD(IZ)**5)**2/(4. *G1(IZ)**1.5) 

C U(IZ)=SQRT(ABS(G1(IZ)-G2(IZ)+G3(IZ))) 

C ENDIF 

157 U(IZ) = SQRT(60. *COORD(IZ) + 100. *COORD(IZ)**6) 

158 UOLD(IZ)= U(IZ) 

160 DO 200 JZ = l.NSNP 

170 A(IZ,JZ) = 0. 

175 B(IZ,JZ) = 0. 

176 C(IZ,JZ) = 0. 

200 CONTINUE 

210 CONTINUE 

C ELEMENTAL DO LOOP TO DETERMINE A MATRIX AND F VECTOR 

214 ALPHA=0. 

215 DO 375 IEL=1,NEL 

220 AE(1,1)=1./ELEN 

230 AE(1,2)=(-1. /ELEN) 

240 AE(2,1)=AE(1,2) 

250 AE(2,2)=AE(1,1) 

260 FS1E(1)=30.*ALPHA*ELEN + 10.*ELEN**2 

270 FS1E(2)=30.*ALPHA*ELEN + 20. *ELEN**2 

272 F1=50.*(ALPHA**6)*ELEN 

274 F2=100.*(ALPHA**5)*(ELEN**2) 

276 F3=125.*(ALPHA**4)*(ELEN**3) 

278 F4=100.*(ALPHA**3)*(ELEN**4) 

280 F5=50.*(ALPHA**2)*(ELEN**5) 

281 F6=100.*ALPHA*(ELEN**6)/7. 

282 F7=25.*(ELEN**7)/16. 

287 FS2E(1)=F1 + F2 + F3 + F4 + F5 + F6 + F7 

290 FS2E(2)=F1 + 2.*F2 + 3.*F3 + 4. *F4 + 5.*F5 + 6.*F6 + 7*F7 

300 DO 370 11=1,2 

310 DO 350 JJ=1.2 

320 IN=ICORR(IEL,II) 
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330 JN=IC0RR(IEL,JJ) 

340 A(IN,JN)=A(IN,JN) - AE(II.JJ) 

350 CONTINUE 

360 FS(IN)=FS1E(II) + FS2E(II) + FS(IN) 

370 CONTINUE 

372 ALPHA=ALPHA + ELEN 

375 CONTINUE 

420 RETURN 

430 END 
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C ******* ** **** * * *** **AAAAAAAA ********* * *** ft A A A A A* A ** A A *** ***** **** * 

C * SUBROUTINE NU2QBI * 

C * * 

C * THIS SUBROUTINE PERFORMS THE ITERATIVE SOLUTION PROCESS FOR * 

C * MAIN PROGRAM NU2QA. a *■ ■ A ■ ; . j^rk***** i a ax 

100 SUBROUTINE NU2QBI(IET) 

102 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 

: UEXT(100),UDIF(100),C00RD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

:ICORR(100,2),ITYPE,NEL,NSNP 

104 DIMENSION WKAREA(40600), DIF(IOO), BE(2,2), USTAR(IOO), FT(IOO), 

:FUE(2),FU(100) 

C SELECT METHOD OF DETERMINING USTAR 

105 PRINT*,'SELECT METHOD OF U* DETERMINATION.' 

106 PRINT*,'1: U* = U' 

107 PRINT*,'2: U* = AVERAGE OF LAST TWO COMPUTED VALUES OF U' 

108 PRINT*,'3: U* = WEIGHTED AVERAGE OF LAST TWO COMPUTED VALUES OF U' 

109 READ(6,*) METHU 

110 IF (METHU. EQ. 3) THEN 

111 PRINT*,'CHOOSE WEIGHTING VALUES A AND B' 

112 READ(6,*) AW,BW 

113 ELSE 

114 CONTINUE 

115 ENDIF 

C SELECT METHOD OF DETERMINING UNSTEADY B MATRIX 

116 PRINT*,'SELECT METHOD OF DETERMINING UNSTEADY B MATRIX. ' 

117 PRINT*,'1: MIDPOINT APPROXIMATION FOR U OVER THE ELEMENT. ' 

119 PRINT*,'2: U LINEARIZED OVER THE LENGTH OF THE ELEMENT.’ 

120 READ(6,*) METHBM 

C SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR 

121 PRINT*,'SELECT METHOD OF DETERMINING UNSTEADY FORCE VECTOR. ’ 

122 PRINT*,'1: MIDPOINT APPROXIMATION’ 

123 PRINT*,'2: 1/4 - 3/4 APPROXIMATION' 

124 PRINT*,'3: LINEAR' 

125 READC6,*) METHFU 

C CALL SUBROUTINE SETIME TO BEGIN TIMING ITERATION PROCESS 

127 CALL SETIME 

C BEGIN ITERATION PROCESS 

128 DO 450 ITER=1,200 

C DETERMINE VALUE OF U* AT EACH NODE AND SET VALUE OF UNSTEADY 
C FORCE VECTOR TO ZERO 

129 DO 138 IU=1,NSNP 

130 FU(IU)=0. 

131 IF (METHU. EQ. 1) THEN 
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132 USTAR(IU)=U(IU) 

133 ELSEIF (METHU. EQ. 2) THEN 

134 USTARCIU)=(U(IU)+UOLD(IU))/2. 

135 ELSE 

136 USTARCIU)=(AW*U(IU)+BW*UOLD(IU))/(AW+Btf) 

137 ENDIF 

138 CONTINUE 

C DETERMINE UNSTEADY ELEMENT B MATRIX 

140 DO 210 IELfI.NEL 

145 IF (METHBM.EQ. 1) THEN 

146 BE(1,1)=(ELEN/3. )*(USTARCIEL)+USTAR(IEL+1)) 

147 B£ 1,2)=CELEN/6. )*(USTAR(IEL)+USTARCIEL+1)) 

148 BEC 2,1)=BEC1,2) 

149 BE(2,2)=BEC1,1) 

150 ELSE 

151 BEC1,1)*CELEN/6. )*(3. *USTAR(IEL) + USTARCIEL+1)) 

152 BEC1,2)=CELEN/6. )*(USTARCIEL) + USTARCIEL+1)) 

153 BEC2,1)=BEC1,2) 

154 BEC2,2)=CELEN/6. )*CUSTAR(IEL) + 3. *USTARCIEL+1)) 

156 ENDIF 

C DETERMINE SYSTEM B MATRIX BY DISTRIBUTING ELEMENT B MATRICES 

C ACCORDING TO THE LOCAL TO GLOBAL CORRESPONDENCE 

157 DO 163 11*1,2 

158 DO 162 JJ*1,2 

159 IN=ICORRCIEL,11) 

160 JN=ICORRCIEL,JJ) 

161 BCIN,JN)=BECII,JJ) + B(IN,JN) 

162 CONTINUE 

163 CONTINUE 

C DETERMINE UNSTEADY ELEMENT FORCE VECTOR 

165 IF CMETHFU.EQ. 1) THEN 

166 FUE(1)=CELEN/2. )*CCUSTARCIEL)+USTARCIEL+1))/2.)**2 

167 FUE(2)=FUE(1) 

169 ELSEIF CMETHFU.EQ. 2) THEN 

171 FUECl)=CELEN/2. )*C3.*USTARCIEL)/4. + USTARCIEL+1)/4.)**2 

172 FUEC2)=CELEN/2. )*CUSTARCIEL)/4. + 3. *USTARCIEL+l)/4. )**2 

173 ELSE 

174 FUEC1)=ELEN*CUSTARCIEL)**2/4. 4USTARCIEL)*USTARCIEL+1)/6. 

: + USTAR EL+1)**2/12. ) 

175 FUEC2)=^EN*CUSTARCIEL)**2/12. + USTARCIEL)*USTARCIEL+1)/6. 

: + USTARCIEL+l)**2/4. ) 

176 ENDIF 

C DETERMINE UNSTEADY SYSTEM FORCE VECTOR BY DISTRIBUTING ELEMENTAL 

C FORCE VECTORS ACCORDING TO THE LOCAL TO GLOBAL CORRESPONDENCE 

177 DO 180 11*1,2 

178 IN=ICORRCIEL,II) 

179 FUCIN)=FUECII) + FUCIN) 
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180 CONTINUE 

210 CONTINUE 

C DETERMINE TOTAL SYSTEM MATRIX 

220 DO 240 IP=1,NSNP 

221 DO 232 JP=1,NSNP 

230 C(IP,JP)=A(IP,JP)+B(IP,JP) 

C RESET B MATRIX TO ZERO 

231 B(IP,JP)=0. 

232 CONTINUE 

* C UPDATE VALUE OF U AT THE PREVIOUS ITERATION 

* 235 UOLD(IP)=U(IP) 

C DETERMINE TOTAL SYSTEM FORCE VECTOR 

236 FT(IP)=FS(IP)+FU(IP) 

240 CONTINUE 

C IMPOSE BOUNDARY CONDITIONS 

241 C(l,l)=l. 

242 C(1,2)=0. 

243 FT(1)=ULBC 

244 IF (ITYPE.EQ. 1) THEN 

245 C(NSNP,NSNP-1)=0. 

246 C(NSNP,NSNP)=1. 

247 FT(NSNP)=URBC 

248 ELSE 

249 FT(NSNP)=FT(NSNP)-URBC 

250 ENDIF 

257 M=1 

260 IDGT=3 

270 IQ=100 

C CALL SUBROUTINE LEQT2F TO SOLVE SET OF LINEAR ALGEBRAIC EQUATIONS 

280 CALL LEQT2F(C,M,NSNP,IQ,FT,IDGT,WKAREA,IER) 

290 DO 310 NEW=1,NSNP 

4 300 U(NEW)=FT(NEW) 

C WRITE(*,*) 'UNEW=',U(NEW) 

* C TEST FOR CONVERGENCE 

306 DIF(NEW)=ABS(U(NEW)-UOLD(NEW)) 

310 CONTINUE 

320 DIFMAX=DIF(1) 

325 NMAX=1 

330 DO 390 IJ=1,NEL 

340 IF (DIF(IJ+1). GE. DIF(IJ)) THEN 

350 DIFMAX=DIF(IJ+1) 

355 NMAX=IJ+1 
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360 ELSE 

370 CONTINUE 

380 ENDIF 

390 CONTINUE 

405 IF (ABS(DIFMAX/U(NMAX)).LT.CONV) THEN 

410 GO TO 460 

420 ELSE 

430 CONTINUE 

440 ENDIF 

450 CONTINUE 

460 CONTINUE 

C CALL SUBROUTINE GETIME TO OBTAIN CPU TIME FOR ITERATION PROCESS 

461 CALL GETIME(IET) 

C OUTPUT HEADER INFORMATION 

462 WRITE(6,464) 

463 WRITE(30,464) 

464 FORMATCIX,'EQUATION: U" + U**2 = 60X + 100X**6') 

465 IF (ITYPE.EQ. 1) THEN 

466 WRITE(6,468) COORDC1),ULBC.COORD(NSNP),URBC 

467 WRITE(30,468) COORD(1),ULBC.COORD(NSNP),URBC 

468 FORMATCIX,'B. C. : U(',F2.0,')=',F2.0,'; U(',F2. 0,')=',F5.0,/) 

469 ELSE 

470 WRITE(6,472) COORD(1),ULBC.COORD(NSNP),URBC 

471 WRITE(30,472) COORDC1),ULBC,COORDCNSNP),URBC 

472 FORMATCIX,'B.C.: UC\F2.0,’)=',F2.0,';DU/DXC',F2. 0,')=',F4.0,/) 

473 ENDIF 

475 IF CMETHU.EQ. 1) THEN 

476 WRITEC6.478) 

477 WRITEC 30,478) 

478 FORMATCIX,'ITERATION METHOD: U*=U’,/) 

479 ELSEIF CMETHU.EQ. 2) THEN 

480 WRITE(6,482) 

481 WRITEC30,482) 

482 FORMATCIX,'ITERATION METHOD: U*=CU+U0LD)/2’,/) 

483 ELSE 

484 WRITEC6,486)AW,BW,AW,BW 

485 WRITEC30,486)AW,BW,AW,BW 

486 FORMATCIX, 1 ITERATION METHOD: U*=C',F4. 1,'*U +',F4. 1,'*UOLD)/C', 
:F4.1 t , +',F4.1, ’)',/) 

487 ENDIF 

488 IF CMETHBM.EQ. 1) THEN 

489 WRITEC6,491) 

490 WRITEC30,491) 

491 FORMATCIX,'METHOD OF B MATRIX INTEGRAL EVALUATION: MIDPOINT',/) 

496 ELSE 

497 WRITEC6.499) 

498 WRITEC30,499) 

499 FORMATCIX,'METHOD OF B MATRIX INTEGRAL EVALUATION: LINEAR’,/) 

500 ENDIF 

501 IF CMETHFU.EQ. 1) THEN 

502 WRITEC6,504) 

503 WRITEC30,504) 
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504 FORMATCIX,'METHOD OF EXCITATION INTEGRAL EVALUATION: MIDPOINT’,/) 

505 ELSEIF (METHFU. EQ. 2) THE 

506 WRITE(6,508) '<* 

507 WRITE(30,508) 

508 FORMATCIX,'METHOD OF EXCITATION INTEGRAL EVALUATION: 1/4-3/4',/) 

509 ELSE 

510 WRITE(6,512) 

511 WRITE(30,512) 

512 FORMATCIX,'METHOD OF EXCITATION INTEGRAL EVALUATION: LINEAR',/) 

513 ENDIF 

515 IF CITER. GE. 200) THEN 

516 WRITEC6,518) 

517 WRITEC30.518) 

518 FORMATCIX,'CONVERGENCE NOT OBTAINED AFTER 200 ITERATIONS.’) 

519 ELSE IF CABSCU(NMAX)).GT. C10.**20).0R.ABSCUCNSNP-1)).GT. 

:CIO.**20)) THEN 

520 WRITEC 6,522) 

521 WRITEC30.522) 

522 FORMATCIX,'SOLUTION PROCESS DIVERGES. ') 

523 ELSE 

524 WRITEC6.526) ITER.NEL 

525 WRITEC30.526) ITER.NEL 

526 FORMATCIX,'CONVERGENCE OBTAINED AFTER ',13,' ITERATIONS USING ', 
: 13,’ ELEMENTS. ’,/) 

530 ENDIF 
540 RETURN 
350 END 
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SUBROUTINE QLEXTB 


C 

c 

c 

c 

c 

c 


********************** 

* 

* 


* 

* 


* THIS SUBROUTINE COMPUTES THE EXACT SOLUTION, U=10X**3, FOR 

* MAIN PROGRAM NU2QB AT THE SPECIFIED NODAL POINTS. 


100 SUBROUTINE QLEXTB 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),U0LD(100), 

:UEXT(100),UDIF(100),C00RD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

: ICORR(100,2),ITYPE,NEL,NSNP 
130 DO 150 NN = l.NSNP 

140 UEXT(NN) = 10. *C00RD(NN)**3 

150 CONTINUE 

160 RETURN * 

170 END 


J 

I 


V 


A 
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C * SUBROUTINE QLOTPT * 

C * THIS SUBROUTINE COMPUTES THE PER CENT ERROR BETWEEN THE EXACT * 

C * AND FEM SOLUTIONS, CPU* FOR THE ITERATION PROCESS, AND PRINTS * 

C * OUT ALL DATA IN TABULAR FORM FOR M AIN PR OGRAM^NU 2 9 A^AND NU 29 B. * 

100 SUBROUTINE QLOTPT(CPUSTAR,IET) 

110 COMMON A(100,100),FS(100),B(100,100),C(100,100),U(100),UOLD(100), 
:UEXT(100),UDIF(100),COORD(100),ELEN,CONV,ELTIME,ULBC,URBC, 

:ICORR(100,2),ITYPE,NEL,NSNP 
115 SUMDIF=0. 

C CALCULATE PER CENT ERROR AT EACH NODE AND SUM THE ABSOLUTE VALUE 
C OF ALL THE ERRORS 

120 DO 150 IK=2,NSNP 

130 UDIF(IK)=100.*(U(IK)-UEXT(IK))/UEXT(IK) 

140 SUMDIF=SUMDIF + ABS(UDIF(IK)) 

150 CONTINUE 

160 UDIF(1)=U(1)-UEXT(1) 

C COMPUTE THE ELAPSED TIME OF THE ITERATION PROCESS 

164 ELTIME=IET*. 000026 

165 WRITE(6,169) ELTIME 

166 WRITE(30,169)ELTIME 

169 FORMAT(IX,'ELAPSED TIME FOR THE ITERATION PROCESS IS ’,F9.4, 

: ' SECONDS.') 

C OUTPUT DATA IN TABULAR FORMAT 

170 WRITE(6,180) 

175 WRITE(30,180) 

180 FORMATC/,IX, X-COORD',4X,'U EXACT',7X,'U FEM',4X,'% DIFF') 

190 WRITE(6,200) (COORD(NP),UEXT(NP),U(NP),UDIF(NP), NP=1,NSNP) 

195 WRITE(30,200) (COORD(NP),UEXT(NP),U(NP),UDIF(NP), NP=1,NSNP) 

200 FORMATC/,2X,F5. 3,4X,F9. 4,3X,F10. 4,4X,F5. 1) 

C CALCULATE CPU* FOR THE ITERATION PROCESS 

205 CPUSTAR=ELTIME*SUMDIF/NSNP 
210 WRITE(6,220) CPUSTAR 

215 WRITE(30,220) CPUSTAR 

220 FORMATC/,IX, CPU* FOR THE ITERATION PROCESS IS ',F9.4,' SECONDS.’) 
230 RETURN 

240 END 
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